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Notes for Students of ESS 223, 
Earth System Physics: 

This monograph on Radiative Transfer provides some core and some supplementary reading ma- 
terial for ESS 223. We will discuss much of the material in the first twenty pages, and the figures 
at the end. The Index beginning on page 168 is also helpful. 

Notes for Students of ESS 236, 
Radiative Transfer and Remote Sensing: 

Yada yada yada. 
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1 Introduction 

This document describes mathematical and computational considerations pertaining to radiative 
transfer processes and radiative transfer models of the Earth system. Our approach is to present 
a detailed derivation of the tools of radiative transfer needed to predict the radiative quantities 
(irradiance, mean intensity, and heating rates) which drive climate. In so doing we begin with 
discussion of the intensity field which is the quantity most often measured by satellite remote 
sensing instruments. Our approach owes much to Bohren and Huffman [1983] (particle scattering), 
Goody and Yung [1989] (band models), and Thomas and Stamnes [1999] (nomenclature, discrete 
ordinate methods, general approach). The nomenclature follows these authors where possible. 
These sections will evolve and differentiate from their original sources as the manuscript takes on 
the flavor of the researchers who contribute to it. 

1.1 Planetary Radiative Equilibrium 

The important role that radiation plays in the climate system is perhaps best illustrated by a simple 
example showing that without atmospheric radiative feedbacks (especially, ironically, the green- 
house effect), our planet's mean temperature would be well below the freezing point of water. 
Earth is surrounded by the near vacuum of space so the only way to transport energy to or from the 
planet is via radiative processes. If E is the thermal energy of the planet, and F^^^ and F^^^ are 
the absorbed solar radiation and emitted longwave radiation, respectively, then 



On timescales longer than about a year the Earth as a whole is thought to be in planetary radiative 
equilibrium. That, is, the global annual mean planetary temperature is nearly constant because the 



dE 
~dt 



pASR _ pOLR 



(1) 
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absorbed solar energy is exactly compensated by thermal radiation lost to space over the course of 
a year. Thus 



i^ASR T7OLR 



(2) 



The total amount of solar energy available for the Earth to absorb is the incoming solar flux (or 
irradiance) at the top of Earth's atmosphere, Fq (aka the solar constant), times the intercepting area 
of Earth's disk which is TTr^. Since Earth rotates, the total mean incident flux kt'^Fq is actually 
distributed over the entire surface area of the Earth. The surface area of a sphere is four times 
its cross-sectional area so the mean incident flux per unit surface area is Fq/A. The fraction of 
incident solar flux which is reflected back to space, and thus unable to heat the planet, is called 
the aplanetary albedo or spherical albedo, R. Satellite observations show that R ^ 0.3. Thus only 
(1 — i?) of the mean incident solar flux contributes to warming the planet and we have 



?ASR 



R)Fe/4: (3) 



Earth does not cool to space as a perfect blackbody (41a) of a single temperature and emissivity. 
Nevertheless the spectrum of thermal radiation F*-""^ which escapes to space and thus cools Earth 
does resemble blackbody emission with a characteristic temperature. The ejfective temperature 
Te of an object is the temperature of the blackbody which would produce the same irradiance. 
Inverting the Stefan-Boltzmann Law (73) yields 

Te = {F^^'^/af^ (4) 

For a perfect blackbody, T = Te. For a planet, the difference between Te and the mean surface 
temperature Tj is due to the radiative effects of the overlying atmosphere. The insulating behavior 
of the atmosphere is more commonly known as the greenhouse effect. 
Substituting (3) and (4) into (2) 

{1-R)FJ4 = aT^ (5) 
il-R)FA''^ (6) 



4a 

For Earth, R ^ 0.3 and Fq ^ 1367 Wm'^. Using these values in (6) yields Te = 255 K. 
Observations show the mean surface temperature T, = 288 K. 



1.2 Fundamentals 

The fundamental quantity describing the electromagnetic spectrum is frequency, v. Frequency 
measures the oscillatory speed of a system, counting the number of oscillations (waves) per unit 
time. Usually v is expressed in cycles-per- second, or Hertz. Units of Hertz may be abbreviated 
Hz, hz, cps, or, as we prefer, s~^. Frequency is intrinsic to the oscillator and does not depend on 
the medium in which the waves are travelling. The energy carried by a photon is proportional to 
its frequency 



E = \iv 



(V) 



1 INTRODUCTION 



3 



where h is Planck's constant. Regrettably, almost no radiative transfer literature expresses quanti- 
ties in frequency. 

A related quantity, the angular frequency to measures the rate of change of wave phase in radi- 
ans per second. Wave phase proceeds through 2n radians in a complete cycle. Thus the frequency 
and angular frequency are simply related 

u = 2niy (8) 

Since radians are considered dimensionless, the units of cu are s^^. However, angular frequency is 
also rarely used in radiative transfer. Thus some authors use the symbol u to denote the element 
of solid angle, as in dco. The reader should be careful not to misconstrue the two meanings. In this 
text we use cu only infrequently. 

Most radiative transfer literature use wavelength or wavenumber. Wavelength, A (m), mea- 
sures the distance between two adjacent peaks or troughs in the wavefield. The universal relation 
between wavelength and frequency is 

\v = c (9) 

where c is the speed of light. Since c depends on the medium, A also depends on the medium. 
The wavenumber z/m~^, is exactly the inverse of wavelength 

1 z/ 

z> = - = - (10) 
A c 

Thus u measures the number of oscillations per unit distance, i.e., the number of wavecrests per 
meter. Using (9) in (10) we find u = u/c so wavenumber D is indeed proportional to frequency 
(and thus to energy). Historically spectroscopists have favored D rather than A or v. Because of this 
history, it is much more common in the literature to find i) expressed in CGS units of cm^ than 
in SI units of m^^. The CGS wavenumber is used analogously to frequency and to wavelength, 
i.e., to identify spectral regions. The energy of radiative transitions are commonly expressed in 
CGS wavenumber units. The relation between i) expressed in CGS wavenumber units (cm^^) and 
energy in SI units (J) is obtained by using (10) in (7) 

E = lOOhci/ (11) 

There is another, distinct quantity also called wavenumber. This secondary usage of wavenum- 
ber in this text is the traditional measure of spatial wave propagation and is denoted by k. 

k = 2ni> (12) 

The wavenumber k is set in Roman typeface as an additional distinction between it and other 
symbols K 

Table 2 summarizes the relationships between the fundamental parameters which describe 
wave-like phenomena. 



'The script k is already used for Boltzmann's constant, absorption coefficients, and vibrational modes 
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Table 2: Wave Parameter Conversion Table' 



Variable 
Units 
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m 


i/ 
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s 
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100c 


271 


271 z/ 
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c 




1 


27CC 


27t 
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lOOA 
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2007rz> 
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100i> 




2007rc 


100ci> 






27rc 






UJ 


271 
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UJ 


2007rc 
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U) 
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kc 


271 
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ck 




2% 




271 


¥ 


20071 






r 


1 


cr 


1 


271 


271 
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lOOcr 


r 


cr 





"The speed of light is c m s . 

^Table entries express the column in terms of the row. 

2 Radiative Transfer Equation 

2.1 Definitions 
2.1.1 Intensity 

The fundamental quantity defining the radiation field is the specific intensity of radiation. Spe- 
cific intensity, also known as radiance, measures the flux of radiant energy transported in a given 
direction per unit cross sectional area orthogonal to the beam per unit time per unit solid angle 
per unit frequency (or wavelength, or wavenumber). The units of Ix are Joule meter^^ second^ ^ 
steradian^ meter^^. In SI dimensional notation, the units condense to Jm~^ s^ sr^ m~^. The SI 
unit of power (1 Watt = 1 Joule per second) is preferred, leading to units of Wm^^ sr~^ m~^ Often 
the specific intensity is expressed in terms of spectral frequency 1^ with units Wm^^ sr~^ Hz~^ or 
spectral wavenumber (also /p) with units Wm^^ sr^^ (cm^^)^^. 

Consider light travelling in the direction Cl through the point r. Construct an infinitesimal 
element of surface area dS intersecting r and orthogonal to Cl. The radiant energy dE crossing dS 
in time dt in the solid angle dVl in the frequency range [u, v + dv\ is related to I^iv, 17) by 

dE = h{v,tl,t,v)dSdtd^dv (13) 

It is not convenient to measure the radiant flux across surface orthogonal to 17, as in (13), when 
we consider properties of radiation fields with preferred directions. If instead, we measure the 
intensity orthogonal to an arbitrarily oriented surface element dA with surface normal n, then we 
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must alter (13) to account for projection of dS onto dA. If the angle between n and Cl is 6 then 



cos 6 = h ■ Cl 



(14) 



and the projection of dS onto dA yields 



dA = cos^d^ 



(15) 



so that 



dE = I^{r, n, t, v) cos e dA dt dVL dv 



(16) 



The conceptual advantage that (16) has over (13) is that (16) builds in the geometric factor required 
to convert to any preferred coordinate system defined by dA and its normal n. In practice dA is 
often chosen to be the local horizon. 

The radiation field is a seven-dimensional quantity, depending upon three coordinates in space, 
one in time, two in angle, and one in frequency. We shall usually indicate the dependence of spec- 
tral radiance and irradiance on frequency by using as a subscript, as in 1^, in favor of the more 
explicit, but lengthier, notation Three of the dimensions are superfluous to climate models 
and will be discarded: The time dependence of 1^ is a function of the atmospheric state and so- 
lar zenith angle and will only be discussed further in those terms, so we shall drop the explicitly 
dependence on t. We reduce the number of spatial dimensions from three to one by assuming a 
stratified atmosphere which is horizontally homogeneous and in which physical quantities may 
vary only in the vertical dimension z. Thus we replace r by z. This approximation is also known 
as a plane-parallel atmosphere, and comes with at least two important caveats: The first is the 
neglect of horizontal photon transport which can be important in inhomogeneous cloud and sur- 
face environments. The second is the neglect of path length effects at large solar zenith angles 
which can dramatically affect the mean intensity of the radiation field, and thus the atmospheric 
photochemistry. 

With these assumptions, the intensity is a function only of vertical position and of direction, 
Iu{z, tt). Often the optical depth r (defined below), which increases monotonically with z, is used 
for the vertical coordinate instead of z. The angular direction of the radiation is specified in terms 
of the polar angle 6 and the azithumal angle (p. The polar angle 6 is the angle between Q, and the 
normal surface n that defines the coordinate system. The specific intensity of radiation traveling 
at polar angle and azimuthal angle (p at optical depth level r in a plane parallel atmosphere is 
denoted by /i/(r, 6, (p). Specific intensity is also referred to as intensity. 

Further simplification of the intensity field is possible if it meets certain criteria. If is not a 
function of position (r), then the field is homogeneous. If is not a function of direction (Cl), then 
the field is isotropic. 

2.1.2 Mean Intensity 

The mean intensity is an integrated measure of the radiation field at any point r. Mean intensity 
is defined as the mean value of the intensity field integrated over all angles. 



u 



(17) 
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The solid angle subtended by Q is the ratio of the area A enclosed by f2 on a spherical surface to the 
square of the radius of the sphere. Since the area of a sphere is 47rr^, there must be An steradians 
in a sphere. It is straightforward to demonstrate that the differential element of area in spherical 
polar coordinates is sin 9 d9 dep. Thus the element of solid angle is 

Q = A/r^ 

dQ = r-^dA 

— r'^ sin 9 d9 d(f) 

= sin^d^d0 (18) 

The field of view of an instrument, e.g., a telescope, is most naturally measured by a solid angle. 

The definition of 1^ (17) demands the radiation field be integrated over all angles, i.e., over all 
471 steradians. Evaluating the denominator demonstrates the properties of angular integrals. The 
denominator of (17) is 

/9=K i'(f>=2n 
/ sin^d^d0 
=0 Jij)=0 

=71 

sin 9 d9 

'0=0 

/0=jt 
sin 9 d9 
=0 

= 27r[-cose]o 

= 2;r[-(-l) - (-1)] 

= 471 (19) 

As expected, there are 471 steradians in a sphere, and 2k steradians in a hemisphere. 

It is convenient to return briefly to the definition of isotropic radiation. Isotropic radiation is, 
by definition, equal intensity in all directions so that the total emitted radiation is simply 471 times 
the intensity of emission in any direction. 

Applying (19) to (17) yields 



An 



h dn (20) 



Ii, has units of radiance, Wm ^ sr ^ Hz ^. If the radiation field is azimuthally independent (i.e., 
ly does not depend on 0), then 



1 r 







2 

Let us simplify (21) by introducing the change of variables 



J^sin^d^ (21) 



u = cos 9 (22) 
du = -sm9d9 (23) 
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This maps 9 E [0, n] into m G [1, —1] so that (21) becomes 



1 



L, du 



I,, du 



(24) 



The hemispheric intensities or half-range intensities are simply the up- and downwelhng com- 
ponents of which the full intensity is composed 



< < 7l/2 
n/2 <6 <Ti 

< M < 1 

-1 < ?i < 



(25a) 
(25b) 



/+(r, /i, = /,(r, (j)) = /,(r, < < 7r/2, (/.) = /,(r, < u < 1, (/.) 
I;{T,fi,(j)) = I^{T,+fi,(f)) = 4(r,7r/2 < < 71,0) = 4(r,-l < M < 0,, 



(26a) 
(26b) 



2.1.3 Irradiance 

The spectral irradiance Fj, measures the radiant energy flux transported through a given surface 
per unit area per unit time per unit wavelength. Although it is somewhat ambiguous, "flux" is used 
a synonym for irradiance, and has become deeply embedded in the literature [Madronich, 1987]. 
Consider a surface orthogonal to the Cl' direction. All radiant energy travelling parallel to Cl' 
crosses this surface and thus contributes to the irradiance with 100% efficiency. Energy travelling 
orthogonal to Q' (and thus parallel to the surface), however, never crosses the surface and does 
not contribute to the irradiance. In general, the intensity Iu{^) projects onto the surface with an 
efficiency cos Q = Cl ■ Cl', thus 



/,, cos 9 dfl 



=271 



I,, cos 6 sin 6 d6 1 



(27) 



=0 J(l)=0 



sin 6 d9. This 



In a plane-parallel medium, this defines the net specific irradiance passing through a given vertical 
level. Note the similarity between (20) and (27). The former contains the zeroth moment of the 
intensity with respect to the cosine of the polar angle, the latter contains the first moment. Also note 
that (27) integrates the cosine-weighted radiance over all angles. If 1^ is isotropic, i.e., Ii, = 1°, 
then Fj, = due to the symmetry of cos 9. 

Let us simplify (27) by introducing the change of variables u = cos 9, du = 
maps 9 e [0, n] into u e [1, -1]: 

pu=—l j'4>='2n 
Fu = / /.M(-dM)d0 

Ju=i J4>=o 

r"M=l l'4>='2% 

I luU du dcj) 

'u=-l J 6=0 



(28) 
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The irradiance per unit frequency, F„, is simply related to the irradiance per unit wavelength, 
Fa. The total irradiance over any given frequency range, [u, v + dz/], say, is Fy ^v. The irradiance 
over the same physical range when expressed in wavelength, [A, A— dA] , say, is Fa dA. The negative 
sign is introduced since — dA increases in the same direction as +dz/. Equating the total irradiance 
over the same region of frequency/wavelength, we obtain 

F^dz/ = -FAdA 

- -^4(1 



-F: 



A 



c A2 

F. = -Fa = -Fa (29) 

Z/"^ c 

Fx = ^F,= -F, (30) 



A 



Thus Ft, and Fa are always of the same sign. 



2.1.4 Actinic Flux 

A quantity of great importance in photochemistry is the total convergence of radiation at a point. 
This quantity, called the actinic flux, F"', determines the availability of photons for photochemical 
reactions. By definition, the intensity passing through a point P in the direction Cl within the solid 
angle d^l is 1^ dVl. We have not multiplied by cos 9 since we are interested in the energy passing 
along Q (i.e., 6 = 0). The energy from all directions passing through P is thus 

F/ = / hdn 

Jin 

= 471 1 (31) 

Thus the actinic flux is simply An times the mean intensity Ii, (20). F;J has units of Wm~^ Hz~^ 
which are identical to the units of irradiance F^ (27). Although the nomenclature "actinic flux" is 
somewhat appropriate, it is also somewhat ambiguous. The "flux" measured by F/ at a point P 
is the energy convergence (per unit time, frequency, and area) through the surface of the sphere 
containing P. This differs from the "flux" measured by F^, which is the net energy transport 
(per unit time, frequency, and area) through a defined horizontal surface. Thus it is safest to use 
the terms "actinic radiation field" for F/ and "irradiance" for F^. Unfortunately the literature is 
permeated with the ambiguous terms "actinic flux" and "flux", respectively. 

The usefulness of actinic flux F^ becomes apparent only in conjunction with additional, species- 
dependent data describing the probability of photon absorption, or photo-absorption. Photo- 
absorption is the process of molecules absorbing photons. Each absorption removes energy (a 
photon) from the actinic radiation field. The amount of photo-absorption per unit volume is pro- 
portional to the number concentration of the absorbing species A^A[m '^l^ the actinic radiation 
field F/, and the efficiency with with each molecule absorbs photons. This efficiency is called the 
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absorption cross-section, molecular cross section, or simply cross-section. The absorption cross- 
section is denoted by a and has units of [m^^]. In the literature, however, values of a usually 
appear in CGS units [cm^^]. To make explicit the frequency-dependence of a we write a{v). If a 
depends significantly on temperature, too (as is true for ozone), we must consider T). 

The probability, per unit time, per unit frequency that a single molecule of species A will absorb 
a photon with frequency in [z/, u + du] is proportional to (z/)a(z/)^. Thus ol{v) is the effective 
cross-sectional area of a molecule for absorption. The absorption cross-section is the ratio between 
the number of photons (or total energy) absorbed by a molecule to the number (or total energy) per 
unit area convergent on the molecule. Let [Wm^'^] be the energy absorbed per unit time, per 
unit frequency, per unit volume of air. Then 

F:{u) = NAF;J{u)a{u) (32) 

where A^a is the number concentration of A. 

Photochemists are interested in the probability of absorbed radiation severing molecular bonds, 
and thus decomposing species AB into constituent species A and B. Notationally this process may 
be written in any of the equivalent forms 

AB + hu — > A + B 
AB + hu A + B 

AB + hu A + B (33) 

Both forms indicate that the efficiency with which reaction (33) proceeds is a function of photon 
energy hu. The second form makes explicit that the photodissociation reaction does not proceed 
unless v < uq, where uq is the photolysis cutoff frequency. In any case, photon energy is conven- 
tionally written hu, rather than the less convenient he/ A. 

The probability that a photon absorbed by AB will result in the photodissociation of AB, and 
the completion of (33), is called the quantum yield or quantum efficiency and is represented by </>. 
The probability (p is dimensionless^. In addition to its dependence on u, (p depends on tempera- 
ture T for some important atmospheric reactions (such as ozone photolysis). We explicitly annotate 
the T-dependence of (p only for pertinent reactions. Measurement of 0(z/) for all conditions and 
reactions of atmospheric interest is an ongoing and important laboratory task. 

The specific photolysis rate coefficient for the photodissociation of a species A is the number 
of photodissociations of A occuring per unit time, per unit volume of air, per unit frequency, per 
molecule of A. In accord with convention we denote the specific photolysis rate coefficient by .J^. 
The units of J,^ are s^^ Hz~^. 

^ Fi{v)a{v)4>{v) 

4%Wa{u)<j){u) ^^^^ 



^The proportionality constant is (hi^)~^, i.e., when the actinic radiation field F;j{v) is converted from energy 
(Jm~^ s~^ Hz^^) to photons (#m^'^ s^^ Hz^^) then F^{v)a{v) is the probability of absorption per second per unit 
frequency 

■'Azimuthal angle and quantum yield are both dimensionless quantities denoted by 0. The meaning of (p should be 
clear from the context. 
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The photon energy in the denominator converts the energy per unit area in to units of photons 
per unit area. The factor of a turns this photon flux into a photo-absorption rate per unit area. 
The final factor, (p, converts the photo-absorption rate into a photodissociation rate coefficient. 
Note that each factor in the numerator of (34) requires detailed spectral knowledge, either of the 
radiation field or of the photochemical behavior of the molecule in question. This complexity is a 
hallmark of atmospheric photochemistry. 

The total photolysis rate coefficient J is obtained by integrating (34) over all frequencies that 
may contribute to photodissociation 



F/(z/)«(z/)0(z/) 



u>u, hi/ 



As mentioned above, evaluation of (35) requires essentially a complete knowledge of the radiative 
and photochemical properties of the environment and species of interest. 

J is notoriously sensitive to uncertainty in the input quantities. Integration errors due to the 
discretization of (35) are quite common. To compute J with high accuracy, regular grids must 
have resolution of ~1 nm in the ultraviolet, [Madronich, 1989]. Much of the difficulty is due to the 
steep but opposite gradients of and that occur in the ultraviolet. High frequency features in a 
worsen this problem for some molecules. 

The utility of J has motivated researchers to overcome these computational difficulties by brute 
force techniques and by clever parameterizations and numerical techniques [Chang et al, 1987; 
Dahlback and Stamnes , 1991; Toon et al, 1989; Stamnes and Tsay, 1990; Petropavlovskikh, 1995; 
Landgraf and Crutzen, 1998; Wild et al, 2000]. It is common to refer to photolysis rate coeffients 
as "J-rates", and to affix the name of the molecule to specify which individual reaction is pertinent. 
Another description for J is the first order rate coefficient in photochemical reactions. For example, 
Jno2 is the first order rate coefficient for 

NOs + hi/ ^<iH^°"^ NO + (36) 

If [NO2] denotes the number concentration of NO2 in a closed system where photolysis is the only 
sink of NO2, then 

= - JN0JNO2] + 5no. (37) 

where 5'no2 represents all sources of NO2. The terms in (37) all have dimensions of m^'^ s^^. The 
first term on the RHS is the photolysis rate of NO2 in the system. 

Figure 1 shows the spectral distribution of actinic flux in a clear mid-latitude summer atmos- 
phere, and the absorption cross-section and quantum yield of NO2. 

Figure 2 shows the vertical distribution of Jno2 [s^^] for the conditions shown in Figure (1). 

In this section we have assumed the quantities F/, a, and are somehow known and therefore 
available to use to compute J. Typically, a and are considered known quantities since they usu- 
ally do not vary with time or space. Models may store their values in lookup tables or precompute 
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Figure 1: (a) Spectral distribution of actinic flux F"^ [#m^^ s^^ |-ini^^] at TOA and at the surface 
for a mid-latitude summer (MLS) atmosphere with a unit optical depth of dust or sulfate in the 
lowest kilometer, (b) Absorption cross section of NO2, aN02 molecule"^], (c) Quantum yield 
ofNO2,0No. (36). 
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Figure 2: Vertical distribution of Jno2 [s^^] (36) for the conditions shown in Figure (1). (a) Ab- 
solute rates, (b) Rates normalized by clean sky rates. 
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Table 3: Actinic Flux Enhancement by Scattering" 



Description p^iu) U{n) F+ /, 



F^ 

CoUimated, F® F®(5(ri - fl©) F® 

non-reflecting 

Isotropic, F® = ° 2^ 

non-reflecting 

CoUimated, 1 ^ F, — 3F, 

reflecting ' 

Isotropic, reflecting F® 1 F®/7r F® — 4F® 

71 

"Terminology: is downwelling irradiance of source, ph{i^) is Lambertian reflectance of surface (or cloud), 
Jy(fi) is resulting intensity field, F+ is upwelling irradiance, is mean intensity, F^J is actinic flux. 

their contributions to (35). The essence of forward radiative problems is to determine Ip so that 
quantities such as Jj, and F" may be determined. In inverse radiative transfer problems, which 
are encountered in much of remote sensing, both J and the species concentration are initially un- 
known and must be determined. We shall continue describing the methods of forward radiative 
transfer until we have tools at our disposal to solve for J^. At that point we shall re- visit the 
inverse problem. 



2.1.5 Actinic Flux Enhancement 

The actinic flux Fj (31) is sensitive to the angular distribution of radiance I,,. Nearly all scattering 
processes diffuse the radiation field, i.e., convert coUimated photons to more isotropic photons. 
Such diffusion causes actinic flux enhancement. It is instructive to examine how the relationship 
between downwelling flux and actinic flux changes in the presence of scattering. Four limiting 
cases may be identified and are summarized in Table 3. The scenarios differ in the isotropicity of 
the downwelling radiance (coUimated or isotropic) and the reflectance Pl('^) (0 or 1) of the lower 
boundary, taken to be a Lambertian surface. All scenarios are driven by the same downwelling 
irradiance, taken to be the direct solar beam. The scenarios are arranged in order of increasing 
actinic flux F^ , shown in the final column. F;j increases with both the number and the brightness 
of reflecting surfaces. The reflectivity of natural surfaces is no more than 90% nor less than 5%**, 
so that the F^/ in Table 3 represent bounds on realistic systems. 

The ColUmated-Non-Reflecting scenario assumes all light travels unidirectionally in a tightly 
coUimated direct beam with irradiance F^, = F®. In this case the radiation field is the delta- 
function in the direction of the solar beam = F®5(i7 — CIq). The closest approximation to 
this scenario in the natural environment is the pristine atmosphere high above the ocean in daylight. 
The actinic flux F^/ = F® follows directly from (31). We define the actinic flux enhancement S 
of a medium as the ratio of the actinic flux to the actinic flux of a coUimated beam with the same 



"^Glaciers are the most reflective surfaces in nature, with Phii') <^ 0.9. Maximum cloud reflectance is <; 0.7. Dark 
forests and ocean have ph{i') ^ 0.05. See discussion in §2.3.2 and Table 4. 
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incident irradiance 

S = F,'/F- (38) 

Table 3 shows 5* ranges from one for the direct solar beam to a maximum of four in a completely 
isotropic radiation field. Thus a coUimated beam is the least efficient configuration of radiant 
energy for driving photochemistry. Scattering processes diffuse the radiation field, and, as a result, 
always enhance the photolytic efficiency of a given irradiance. Absorption in a medium (i.e., in 
the atmosphere or by the surface) always reduces and may lead to S < 1. 

The Isotropic -Non-Reflecting scenario assumes isotropic downwelling radiance above a com- 
pletely black surface. Under these conditions, the actinic flux is twice the incident irradiance 
because the photons are evenly distributed over the hemisphere rather than coUimated. Moderately 
thick clouds (r ^ 3) over a dark surface such as the ocean create a radiance field approximately 
hke this. However, because clouds are efficient at diffusing downwelling irradiance, they are also 
efficient reflectors and this significantly reduces the incident irradiance relative to the total 
extraterrestrial irradiance F®. Whether photochemistry is enhanced or diminished beneath real 
clouds depends on whether the actinic flux enhancement factor F;J = 2 compensates the reduced 
sub-cloud insolation due to photons up-scattering off the cloud and back to space. 

The CoUimated-Reflecting scenario is a very important limit in nature because the net effect 
on photochemistry can approach the theoretical photochemical enhancement of S' = 3. In this 
limit, coUimated downwelling radiation and diffuse upwelling radiation combine to drive pho- 
tochemistry from both hemispheres. A clear atmosphere above bright surfaces (clouds, desert, 
snow) approaches this limit. These conditions describe a large fraction of the atmosphere, which 
is 50-60% cloud-covered. Moreover, the incident flux is not attenuated by cloud transmission, so 
F~ ^ F®. The relatively high frequency of occurance of this scenario, combined with the large 
photochemical enhancement of 5 = 3, are unequalled by any other scenario. 

We may also define the actinic flux efficiency F of a medium as the actinic flux relative to the 
actinic flux of an isotropic radiation field with the same incident irradiance 

F ^ ^ (39) 

Clearly < F < 1. Table 3 shows F = 0.25 for a coUimated beam, and F = 1 for isotropic 
radiation. Just as scattering of the solar beam is required to increase F above 0.25, absorption must 
be present to reduce F beneath 0.25. 

The Isotropic-Reflecting scenario shows that an isotropic radiation field is most efficient for 
driving photochemistry. The maximum four-fold increase in efficiency relative to the coUimated 
field arises from the radiation field interacting with the particle from all directions, rather than from 
one direction only. This is geometrically equivalent to multiplying the molecular cross section 
by a factor of four, the ratio between the surface and cross-sectional areas of a sphere. Note 
that photochemistry itself is driven by molecular absorption which reduces F^ from the values in 
Table 3. 

In summary, we have learned that a reactant molecule with a spherically symmetric field of 
influence receives photochemical radiation much like Earth receives solar irradiance. In both cases 
the coUimated beam intercepts one fourth of the total area of matter (molecule or Earth), while an 
equal flux of diffuse (or diurnal average) irradiance impinges on four times as much area. Since, 
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in the geometric limit, absorption probability depends upon area, not direction, coUimated beams 
have one fourth the photochemical potential as isotropic radiation. 



2.1.6 Energy Density 

Another quantity of interest is the density of radiant energy per unit volume of space. We call this 
quantity the energy density Uy. The energy density is the number of photons per unit volume in 
the frequency range [z/, u + dz/] times the energy per photon, hz/. U,, is simply related to the actinic 
flux (31) and thus to the mean intensity 1^. 

= [ dU, 

J An 

An - 

= —h (40) 

c 

The units of Uy are J m~'^ Hz^^. 



2.1.7 Spectral vs. Broadband 

Until now we have considered only spectrally dependent quantities such as the spectral radiance ly, 
spectral irradiance Fy, spectral actinic flux F/, and spectral energy density Uy. These quantities are 
called spectral and are given a subscript of u. A, or u because they are expressed per unit frequency, 
wavelength, or wavenumber, respectively. Each spectral radiant quantity may be integrated over 
a frequency range to obtain the corresponding band-integrated radiant quantity. Band- integrated 
radiant fields are often called narrowband or broadband Depending on the size of the frequency 
range, broadband radiant are obtained by integrating over all frequencies: 



F = I Fju)du 



F^ = / F:!(u)du 



F'' = F^{v)dy 
Jo 

POO 

J = Jyiv^dv 



U = / Uy{u)du 
Jo 

POD 

h = hyMdu 



2.1.8 Thermodynamic Equilibria 

Temperature plays a fundamental role in radiative transfer because T determines the population of 
excited atomic states, which in turn determines the potential for thermal emission. Thermal emis- 
sion occurs as matter at any temperature above absolute zero undergoes quantum state transitions 
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from higher energy to lower energy states. The difference in energy between the higher and lower 
level states is transferred via the electromagnetic field by photons. Thus an important problem in 
radiative transfer is quantifying the contribution to the radiation field from all emissive matter in a 
physical system. For the atmosphere the system of interest includes, e.g., clouds, aerosols, and the 
surface. 

To develop this understanding we must discuss various forms of energetic equilibria in which 
a physical system may reside. Earth (and the other terrestrial planets. Mercury, Venus, and Mars) 
are said to be in planetary radiative equilibrium because, on an annual timescale the solar energy 
absorbed by the Earth system balances the thermal energy emitted to space by Earth. Radiation 
and matter inside a constant temperature enclosure are said to be in thermodynamic equilibrium, 
orTE. 

Radiation in thermodynamic equilibrium with matter plays a fundamental role in radiation 
transfer. Such radiation is most commonly known as blackbody radiation. Kirchoff first deduced 
the properties of blackbody radiation. 

Thermodynamic equilibrium (TE) is an idealized state, but, fortunately, the properties of radi- 
ation in TE can be shown to apply to a less restrictive equilibrium known as local thermodynamic 
equilibrium, or LTE. 



2.1.9 Planck Function 

The Planck function B,, describes the intensity of blackbody radiation as a function of temperature 
and wavelength 

^^(^'-) = ^^(^Wi^ (41a) 
2hc^ 

^^(^'^) = A^(^S^Al^ ^^^^^ 

Blackbodies emit isotropically — this considerably simplifies thermal radiative transfer. The cor- 
rect predictions (41a) resolved one of the great mysteries in experimental physics in the late 19th 
century. In fact, this discovery marked the beginning of the science of quantum mechanics. 

The relations (41a) and (41b) predict slightly different quantities. The former predicts the 
blackbody radiance per unit frequency, while the latter predicts the blackbody radiance per unit 
wavelength. Of course these quantities are related since the blackbody energy within any given 
spectral band must be the same regardless of which formula is used to describe it. Expressed 
mathematically, this constraint means 

Bydv = -BxdX (42) 



Once again, the negative sign arises as a result of the opposite senses of increasing frequency 
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versus increasing wavelength. We may derive (41b) from (41a) by using (9) in (42) 

= —B\ — 
au 

Bxc 

= -B^ = ^B^ (43) 

2 

Bx = -B, = ^B, (44) 
c X'' 

These relations are analogous to (30). 

The Planck function (41) has interesting behavior in both the high and the low energy photon 
limits. In the high energy limit, known as Wien's limit, the photon energy greatly exceeds the 
ambient thermal energy 

hm > kT (45) 

In Wien's limit, (41) becomes 

S.(r,z.)=^e-Wk^ (46a) 
Bx{T,X)~e-'^^/'^^ (46b) 

In the very low energy limit, known as the Rayleigh-Jeans limit, the photon energy is much 
less than the ambient thermal energy 

huu < kT (47) 

Thus in the Rayleigh-Jeans limit the arguments to the exponential in (41) are less than 1 so the 
exponentials may be expanded in Taylor series. Starting from (41a) 



2hi/3 / hz/ ^ ^ 
2hi/3kT 



B,{T,u) ^ _(l + _-l 



c 

Similar manipulation of (41b) may be performed and we obtain 



c^hz/ 



2u\T 

B,{T,u)^^^ (49a) 
9pkT 

5a(T,A)^^ (49b) 
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The frequency of extreme emission is obtained by taking the partial derivative of (41a) with 
respect to frequency with the temperature held constant 



dv 



X ^ X ( ST/^fpWkT _ 1^ _ ^3 JLpWkT^ 



To solve for the frequency of maximum emission, z/m, we set the RHS equal to zero so that one or 
more of the LHS factors must equal zero 



2hz/2^ [3(6^^^^/'^^ - 1) - ii^Meh-M/kTj 







3(^ghi.M/kT _ 1) _ ^^ghi^M/kT ^ g 



kT 

hw 



ghi/M/kT ^ 2^ghj/M/kT _ 1) 

- - 3(1-6-'^'^^*/'^^) (50) 



kT 

An analytic solution to (50) is impossible since z/m cannot be factored out of this transcendental 
equation. If instead we solve x = 3(1 — e^"^) numerically we find that x ~ 2.8215 so that 



M 



2.82 



kT 

Z/M ^ 2.82kT/h 

^ 5.88 X 10^° T Hz (51) 

where the units of the numerical factor are Hz Thus the frequency of peak blackbody emission 
is directly proportional to temperature. This is known as Wien 's Displacement Law. 

A separate relation may be derived for the wavelength of maximum emission Am by an analo- 
gous procedure starting from (41b). The result is 

Am ^ 2897.8/T |im (52) 

where the units of the numerical factor are |J.mK. Note that (51) and (52) do not yield the same 
answer because they measure different quantities. The wavelength of maximum emission per unit 
wavelength, for example, is displaced by a factor of approximately 1 .76 from the wavelength of 
maximum emission per unit frequency. 

It is possible to use (5 1) to estimate the temperature of remotely sensed surfaces. For example, a 
satellite-borne tunable spectral radiometer may measure the emission of a newly discovered planet 
at all wavelengths of interest. Assuming the wavelength of peak measured emission is Am l-im. 
Then a first approximation is that the planetary temperature is close to 2897. 8/ Am- 



2.1.10 Hemispheric Quantities 

In climate studies we are most interested in the irradiance passing upwards or downwards through 
horizontal surfaces, e.g., the ground or certain layers in the atmosphere. These hemispheric irradi- 
ances measure the radiant energy transport in the vertical direction. These hemispheric irradiances 
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depend only on the corresponding hemispheric intensities. Let us assume the intensity field is 
azimuthally independent, i.e., ly = Iu{d) only. Then the azimuthal contribution to (28) is 271 and 

/•u=l 

Fy = 2n lyudu 

Ju=~l 

(/•u=0 /•u=l \ 

/ hudu+ / lyudu] (53) 
Ju=~l Ju=0 J 

We now introduce the change of variables /i = \u\ = \ cos9\. Referring to (23) we find 

I ni f cos9 : < 9 < n/2 . , 

a=\cos9\ = <. „ ,„ „ ' (54a) 

' ' \ -cos9 : 71/2 < 9 <% 

. . ( u : < u < 1 N 
H = \u\ = < - „ (54b) 

Most formal work on radiative transfer is written in terms of fi rather than u or 9. Substituting (54) 
into (53) 

/ /•/i=0 /•/i=l 

Fy = 271 / /^.(-/i) (-d/i) + / lufidfi 
= 271 I / Jy/i d/i + / d/i 

\Jti=l JtJ.=0 

= 271 I — / Jy/i d/i + / d/i 

V J/x=0 ~'At=0 

= -K + K 

= F^-K (55) 
where we have defined the hemispheric fluxes or half -range fluxes 

F+ = 2ti[ Iy{+fi)^dfi (56a) 
Jo 

F- = 211 I /^(-/i)/td/i (56b) 
The hemispheric fluxes are positive-definite, and their difference is the net flux. 



F, = F+ - F- (57) 



F = F^ - F~ 



(58) 



The superscripts and ~ denote upwelling (towards the upper hemisphere) and downwelling (to- 
wards the lower hemisphere) quantities, respectively. The net flux F^, is the difference between 
the upwelling and downwelling hemispheric fluxes, F+ and F^ , which are both positive-definite 
quantities. 
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The hemispheric flux transport in isotropic radiation fields is worth examining in detail since 
this condition is often met in practice. It will be seen that isotropy considerably simplifies many 
of the troublesome integrals encountered. When I,, has no directional dependence (i.e., it is a 
constant) then = Fj7 (58). Thus the net radiative energy transport is zero in an isotropic 
radiation field, such as a cavity filled with blackbody radiation. 

Let us compute the upward transport of radiation B^{T) (41a) emitted by a perfect blackbody 
such as the ocean surface. Since Bj,(T) is isotropic, the intensity may be factored out of the 
definition of the up welling flux (56a) and we obtain 



2% / Byjidjj, 



pi 


= 2kB^ 


/id/i 


Jo 




= 2nB,^ 


fi=0 


= 2'!iB,{\ 


-0) 


= nB^ 





(59) 



Thus the upwelling blackbody irradiance tranports n times the constant intensity of the radiation. 
Given that the upper hemisphere contains 2n steradians, one might naively expect the upwelling 
irradiance to be 2'!iBy. In fact the divergence of blackbody radiance above the emitting surface 
is 2'kB^. But the vertical flux of energy is obtained by cosine- weighting the radiance over the 
hemisphere and this weight introduces the factor of | difference between the naive and the correct 
solutions. 



2.1.11 Stefan-Boltzmann Law 

The frequency-integrated hemispheric irradiance emanating from a blackbody of great interest 
since it describes, e.g., the radiant power of most surfaces on Earth. Although we could integrate 
the By (41a) directly to obtain the broadband intensity B = B^^u) du, it is traditional to 
integrate B^, first over the hemisphere (59). By proceeding in this order, we shall obtain the total 
hemispheric blackbody irradiance in terms fundamental physical constants and the temperature 
of the body. 

POO 

F+ = n Bydv 
Jo 

- c2(eWkT_i) 
27rh r 

,,W^T _ 1 (60) 
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To simplify (60) we make the change of variables 

hu 

^ = kT 
kTx 

" = IT 

du = —— ax 
h 

dx = —— dp 
kT 

This change of variables maps v G [0, oo) to x G [0, oo). Substituting this into (60) we obtain 

2*'^' d. (61) 



c^h'^ 7q — 1 

The definite integral in (61) is n'^/lb. Proving this is a classic problem in mathematical physics 
which involves the Riemann zeta function (and thus prime number theory), the Gamma function, 
and contour integration. The procedure used to obtain this result is interesting so we briefly sum- 
marize it here. The Riemann zeta function ((x) for real x > 1 may be defined as 

C(x) = — / ^— d« (62) 
Comparing (61) with the Riemann zeta function definition (62), we see that x = 4, i.e, 

n = ^r(4)C(4) (63) 

The integral (62) is analytically solvable for the special case of integers x = n. We may 
transform the integrand from a rational fraction into a power series using algebraic manipulation: 



e" - 1 





e 
X - 
e 


e" - 1 


r— 1 - 

u e 


-u 


1 -e- 


■u 




■u ^ 



1 - e-" 



(64) 



The integration limits in (62) are [0, oo) so it is always true that e~" < 1 in (64), and thus in the 
integrand of (62). 

Recall that the sum of an infinite power series with initial term and ratio r is 



> a^iT^ = lim Oo + aor + a^r'^ + ■ ■ ■ + a^r^ ^ + a^r^ 

ao/(l - r) for |r| < 1 (65) 



fc— >oo 
fc=0 
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Hence the last term in (64) is the sum of a power series (65) with initial term oq = 1 and ratio 



1 



1 - e-" 



k=0 



-ku 



Substituting (66) into (64) we obtain 



u 



x-l 



e" - 1 



u 



-ku 



k=0 



k=0 

CO 



-(k+l)u 



-ku 



k=l 



where the last step shifts the initial index to from zero to one. 

Using the inifinite series representation (67) for the integrand of (62) yields 



ax) 



r(x) 



-ku 



k=l 



du 



Integration and addition are commutative operations. Interchanging their order yields 



Fx 



^ oo 



k=l 



^x-l^-ku 



We change variables from u E [0, +oo] to y G [0, +oo] with 

y = ku 

dy = kdu 

u = y/k 

du = k^^ dy 



so that (68) becomes 



T(x 



k=l 

oo 



-^k-^ dy 



-i^5^A;-(-^)fc-i / y^-'e-ydy 

^ oo 

f(^E*-'irWl 

^ ' k=l 

oo 



k=l 



(66) 



(67) 



(68) 



(69) 



(70) 
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where we replaced the integral in brackets with the Gamma function it defines. Hence the Riemann 
zeta function of a positive integer n is the sum of the reciprocals of the postive integers to the 
power n. 

Contour integration in the complex plane gives analytic closed-form solutions to T.fk~"' (70) 
for positive, even integers n. Our immediate concern is n = 4. Consider the complex function 
[Carrier etal, 1983, p. 97] 



TTCOtTTZ 

fiz) = — — (VI) 



This function 



1 . Is analytic throughout the complex plane 

2. Has first order poles at all integer values on the real axis (except the origin) 

3. Has a fifth order pole at the origin 

4. Satisfies lim|/j|_^oo fiz) = where z = Re^^ 

Therefore (71) obeys the residue theorem for suitably chosen contours. In other words, fxm. 
Having shown 

C(4) = 71790 (72) 
Using r(4) = 3! = 3 X 2 X 1 = 6 and C(4) = 7^790 from (72), Equation (63) becomes 

F+ = — — X 6 X — X 
^ c%3 90 



15c%3- 

This is known as the Stefan-Boltzmann Law of radiation, and is usually written as 

F+ = aT^ where (73) 

' " 15^ (^^^ 

a is known as the Stefan-Boltzmann constant and depends only on fundamental physical constants. 
The value of a is 5.67032 x 10^*^ Wm^^ K^''. The thermal emission of matter depends very 
strongly (quartically) on T (73). This rather surprising result has profound implications for Earth's 
chmate. 

We derived (73) directly so that the Stefan-Boltzmann constant would fall naturally from 
the derivation. For completeness we now present the broadband blackbody intensity B 

POO 

B = B^{iy)du 
Jo 



15c2h3- 

= F^/n = (7X^71 (75) 

The factor of n difference between B and F^ is at first confusing. One must remember that B is the 
broadband (spectrally-integrated) intensity and that F^ is the broadband hemispheric irradiance 
(spectrally-and-angularly-integrated). 
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2.1.12 Luminosity 

The total thermal emission of a body (e.g., star or planet) is called its luminosity, L. The luminosity 
is thermal irradiance integrated over the surface-area of a volume containing the body. The lumi- 
nosity of bodies with atmospheres is usually taken to be the total outgoing thermal emission at the 
top of the atmosphere. If the time-mean thermal radiation field of a body is spherically symmetric, 
then its luminosity is easily obtained from a time-mean measurement of the thermal irradiance 
normal to any unit area of the surrounding surface. This technique is used on Earth to determine 
the intrinsic luminosity of stars whose distance Dq is known (e.g., through parallax) including our 
own. 

Le = 47iD|Fq (76) 

The meaning of the solar constant Fq is made clear by (76). Fq is the solar irradiance that would 
be measured normal to the Earth-Sun axis at the top of Earth's atmosphere. The mean Earth- 
Sun distance is ~ 1.5 x 10^^ m and Fq ^ ISeTWm^^ so Lq = 3.9 x lO^^W. Note that L has 
dimensions of power, i.e., J s^^. 

An independent means of estimating the luminosity of a celestial body is to integrate the surface 
thermal irradiance over the surface area of the body. For planets without atmospheres, L is simply 
the integrated surface emission. Assuming the effective temperature (4) of our Sun is Tg, the radius 
Tq at which this emission must originate is defined by combining (73) with (76) 



^2 

'q 

re 

For the Sun-Earth system, Tq ^ 5800 K so = 6.9 x 10*^ m. Solar radiation received by Earth 
appears to originate from a portion of the solar atmosphere known as the photosphere. 



= iTlD^F^ 

DIFq 




(77) 



2.1.13 Extinction and Emission 

Radiation and matter have only two forms of interactions, extinction and emission. Lambert^, 
first proposed that the extinction (i.e., reduction) of radiation traversing an infinitesimal path ds is 
linearly proportional to the incident radiation and the amount of interacting matter along the path 

= —k(u)I„ Extinction only (78) 

ds 

Here k(i') is the extinction coefficient, a measurable property of the medium, and s is the absorber 
path length, k^u) is proportional to the local density of the medium and is positive-definite. 

The term extinction coefficient and the exact definition of k{u) are somewhat ambiguous un- 
til their physical dimensions are specified. Path length, for example, can be measured in terms 



^Beer is usually credited with first formulating this law, but actually the similarly-named Rougher deserves the 
credit. 
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of column mass path M[kgm^^], number path of molecules A^[#m"^], and geometric distance 
s [m]. Each path measure has a commensurate extinction coefficient: the mass extinction coef- 
ficient ipe [m^kg^^] (i.e., optical cross-section per unit mass), the number extinction coefficient 
Ke [m^ molecule"^] (i.e., optical cross-section per molecule), the volume extinction coefficient 
fee [m^ m"^]=[m"^] (i.e., optical cross-section per unit concentration). These coefficients are inter- 
related by 



[m^ molecule^ ^] 



p [kg m-3] 

ipc [m^kg"^] X M [kgmol^^] 



A/^ [molecule mol ^] 
fee [m^^] = Kc [m^molecule^^] X n [molecule m~^] (79) 

We will develop the formalism of radiative transfer in this chapter in terms using the geometric 
path and volume extinction coefficient (s, ke) formalism. Our intent is to be concrete, rather than 
leaving the choice of units unstated. However, there is nothing fundamental about (s, kc). In 
Section 4.1.1 we state our preference for working in mass-path units (M, 'ipc). 

Extinction includes all processes which reduce the radiant intensity. As will be described 
below, these processes include absorption and scattering, both of which remove photons from 
the beam. Similarly the radiative emission is also proportional to the amount of matter along the 
path 

= kMSu Emission only (80) 

as 

where S„ is known as the source function. The source function plays an important role in radiative 
transfer theory. We show in §2.2.1 that if 3,^ is known, then the full radiance field Jj, is deter- 
mined by an integration of with the appropriate boundary conditions. Emission includes all 
processes which increase the radiant intensity. As will be described below, these processes include 
thermal emission and scattering which adds photons to the beam. Determination of k(i'), which 
contains all the information about the electromagnetic properties of the media, is the subject of 
active theoretical, laboratory and field research. 

Extinction and emission are linear processes, and thus additive. Since they are the only two 
processes which alter the intensity of radiation. 



1 dh 



-h + (81) 



kc ds 

Equation (81) is the equation of radiative transfer in its simplest differential form 
2.1.14 Optical Depth 

We define the optical path f between points Pi and P2 as 

Pi 



f{Pi,P2) = [ ' kods 



'Pi 

df = keds (82) 
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The optical path measures the amount of extinction a beam of light experiences traveling between 
two points. When f > 1, the path is said to be optically thick. 

The most frequently used form of optical path is the optical depth. The optical depth r is the 
vertical component of the optical path f, i.e., r measures extinction between vertical levels. For 
historical reasons, the optical depth in planetary atmospheres is defined r = at the top of the 
atmosphere and r = r* at the surface. This convention reflects the astrophysical origin of much 
of radiative transfer theory. Much like pressure, r is a positive-definite coordinate which increases 
monotonically from zero at the top of the atmosphere to its surface value. Consider the optical 
depth between two levels Z2 > zi, and then allow Z2 ^ oo 

t{zi,Z2) = / hdz' 

Jz'=zi 
pz'=oo 

riz, oo) = / kgdz' (83) 

Jz'=z 
z'=oo 

K <^z' (84) 

J=z 

Equation (84) is the integral definition of optical depth. The differential definition of optical depth 
is obtained by differentiating (84) with respect to the lower limit of integration and using the 
fundamental theorem of differential calculus 

dr = [kc{oo) — kc{z)]dz 

= -kg{z)dz (85) 

where the second step uses the convention that fce(oo) = 0. By convention, r is positive-definite, 
but (85) shows that dr may be positive or negative. If this seems confusing, consider the anal- 
ogy with atmospheric pressure: pressure increases monotonically from zero at the top of the at- 
mosphere, and we often express physical concepts such as the temperature lapse rate in terms of 
negative pressure gradients. 

2.1.15 Geometric Derivation of Optical Depth 

The optical depth of a column containing spherical particles may be derived by appealing to in- 
tuitive geometric arguments. Consider a concentration of A^m^^ identical spherical particles of 
radius r residing in a rectangular chamber measuring one meter in the x and y dimensions and of 
arbitrary height. If the chamber is uniformly illuminated by a collimated beam of sunlight from 
one side, how much energy reaches the opposite side? 

For this thought experiment, we will neglect the effects of scattering^'. Moreover, we will 
assume that the particles are partially opaque so that the incident radiation which they do not absorb 
is transmitted without any change in direction. Finally, assume the particles are homogeneously 
distributed in the horizontal so that the radiative flux F{z) is a function only of height in the 
chamber. Let us denote the power per unit area of the incident collimated beam of sunlight as Fq. 
Our goal is to compute F{z) as a function of and of r. 




*The effects of diffraction will not be explicitly included either, but are implicit in the assumption of the horizontal 
homogeneity of the radiative flux. 



2 RADIATIVE TRANSFER EQUATION 



27 



Since each particle absorbs incident energy in proportion to its geometric cross-section, the 
total absorption of radiation per particle is proportional to nr^Qc, where Qe = I for perfectly 
absorbing particles. If Qc < 1, then each particle removes fewer photons than suggested by its 
geometric size^. Conversely, if Qc > 1 each particle removes more photons than suggested by 
its geometric size. Each particle encountered removes nr'^QeF{z) Wm~^ from the incident beam. 
The maximum flux which can be removed from the beam is, of course, Fq. 

In a section of height Az, the coUimated beam passing through the chamber will encounter 
a total of NAz particles. The number of particles encounted times the flux removed per particle 
gives the change in the radiative flux of the beam between the entrance and the rear wall 



AF = -nr^Q^FNAz 

^ = -nr^Q.NAz (86) 
r 

If we take the limit as Az then (86) becomes 

^dF = -nr^Q^Ndz 
F 

d(lnF) = -nr^Q^Ndz (87) 

Let us define the volume extinction coefficient k as 

k = nr^Q^N (88) 

Finally we define the optical depth r in terms of the extinction 

T = kAz (89) 
r = nr^Q^NAz (90) 

The dimensions of k are inverse meters [m^^], or, perhaps more intuitively, square meters of effec- 
tive surface area per cubic meter of air [m^ m^^]. Therefore r is dimensionless. It can be helpful 
to remember that All quantities which compose r are positive by convention, therefore r itself is 
positive-definite. 

We are now prepared to solve (87) for F(t). From the theory of first order differential equa- 
tions, we know that F must be an exponential function whose solution decays from its initial value 
with an e-folding constant of r 

F(r) = F^e-^ (91) 

Thus, in the limit of geometrical optics, the optical depth measures the number of e-foldings un- 
dergone by the radiative flux of a coUimated beam passing through a given medium. This result is 
one form of the extinction law. 

The exact value of Qc{r, A) depends on the composition of the aerosol. However, there is a 
hmiting value of Qe as particles become large compared to the wavelength of light. 

lim Qe = 2 (92) 



^Although intuition suggests a spherical particle should not remove more energy from a coUimated light beam than 
it can geometrically intercept, this is not the case. As discussed later, diffraction around the particle is important. In 
fact, for r » A, Qe — > 2. However, most of the extinction due to diffraction occurs as scattering, not absorption. 



2 RADIATIVE TRANSFER EQUATION 



28 



Thus particles larger than 5-10 |im extinguish twice as much visible light as their geometric cross- 
section suggests. 

To gain more insight into the usefulness of the optical depth, we can express r (90) in terms of 
the aerosol mass M, rather than number concentration N. For a monodisperse aerosol of density 
p, the mass concentration is 



M 



4 



-nr^pN 



(93) 



If we substitute 



into (90) we obtain 



3M 



r 



4rp 
4rp 



(94) 



Typical cloud particles have r ~ 10 |im so that for visible solar radiation with A ~ 0.5 |im we may 
employ (92) to obtain 

SMAz 

Thus r increases linearly with M for a given r. Note however, that a given mass M produces an 
optical depth that is inversely proportional to the radius of the particles! 



2.1.16 Stratified Atmosphere 

We obtain the radiative transfer equation in terms of optical path by substituting (82) into (81) 

^ = -h + S, (96) 
dr 

A stratified atmosphere is one in which all atmospheric properties, e.g., temperature, density, vary 
only in the vertical. As shown in the non-existant Figure, the photon path increment ds at polar 
angle in a stratified atmosphere is related to the vertical path increment dz by dz = cos 9 ds or 
ds = u^^ dz. In other words, the optical path traversed by photons is proportional to the vertical 
path divided by the cosine of the trajectory. 

df = u~^ dr (97a) 
dr = n df (97b) 

Substituting (97a) into (96) we obtain 

= -h + S, (98) 

dr 

This is the differential form of the radiative transfer equation in a plane parallel atmosphere and 
is valid for all angles. The solution of (98) is made difficult because Su depends on I,^. A more 
tractable set of equations may be obtained by considering the form of the boundary conditions. 
For many (most) problems of atmospheric interest, we know I„ over an entire hemisphere at each 
boundary of a "slab". Considering the entire atmosphere as a slab, for example, we know that. 
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at the top of the atmosphere, sunhght is the only incident intensity from the hemisphere contain- 
ing the sun. Or, at the surface, we have constraints on the upwelling intensity due to thermal 
emission or the surface reflectivity. When combined, these two hemispheric boundary conditions 
span a complete range of polar angle, and are thus sufficient to solve (98). However, in practice 
it is difficult to apply half a boundary condition. Moreover, we are often interested in knowing 
the hemispheric flows of radiation because many instruments (e.g., pyranometers) are designed 
to measure hemispheric irradiance and many models (e.g., climate models) require hemispheric 
irradiance to compute surface exchange properties. For these reasons we will decouple (98) into 
its constituent upwelling and downwelling radiation components. 

Using the definitions of the half-range intensities (25) in (98) we obtain 

-u^^ = 1- -S- Q<e < 71/2, n = COS0 > (99a) 
dr 

d/+ 

-u—^ = 1+ -S^ Ti/2 <9 <K,u = cose <0 (99b) 
dr 

where we have simply multiplied (98) by — 1 in order to place the negative sign on the LHS for 
reasons that will be explained shortly. The definitions of and 5*+ are exactly analogous to (25). 

We now change variables from u to n (54). Replacing m by in (99a) is allowed since fi = 
u > in this hemisphere. In the upwelling hemisphere where 7r/2<6'<7r, n<Oso that fi = —u 
(54). This negative sign cancels the negative sign on the LHS of (99b), resulting in 

= C - s; (looa) 

= /+ - S+ (100b) 

These are the equations of radiative transfer in slab geometry for downwelling (0 < 9 < Ti/2) 
and upwelling (71/2 < 9 < n) intensities, respectively. The only mathematical difference between 
(100a) and (100b) is the negative sign. A helpful mnemonic is that the negative sign on the LHS 
is associated with while the implicit unary positive sign is associated with /+. Of course the 
and Si, terms on the RHS are prefixed with opposited signs since they represent opposing, but 
positive definite, physical processes (absorption and emission). 

Equation (99) states that I"^ depends explicitly only on and on 5*^, but has no explicit 
dependence on or on S^. Thus it may appear that /"*" and /~ are completely decoupled from 
eachother. However, we shall see that in problems involving scattering, S*^ depends explicitly on 
P because scattering may change the trajectory of photons from upwelling to downwelling and 
visa versa. By coupling to scattering allows the entire radiance field to affect the radiance 
field at every point and in every direction (modulo the speed of light, of course). Thus scattering 
changes the solutions to (100) from being locally-dependent to depending on the global radiation 
field. 

As a special case of (98), consider a stratified, non-scattering atmosphere in thermodynamic 
equilibrium. Then the source function equals the Planck function Su = By = By{T) and we have 

u'^ = -Iy + By (101) 

dr 

Equation (101) is the basis of radiative transfer in the thermal infrared, where scattering effects are 
often neghgible. The solution to (101) is described in §2.2.2. 
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2.2 Integral Equations 
2.2.1 Formal Solutions 

It is useful to write down the formal solution to (98) before making additional assumptions about 
the form of the source function S^. 

U—— = -U + Ou 
(XT 

u dli, = —Ii, dr + Sy dr 
u diu + dr = Si, dr 

d/^ + — dr = M^^^^dr (102) 
u 

Multiplying (102) by the integrating factor e"^/" 

e^/'^d/^ + ^dr = M-V/"^^dr 

u 

d(e"/"J^) = M-ie"/"5^dr 
d(e^/"4) 



dr 



u 



^e"/"^^ (103) 



The LHS side of (103) is a complete differential. The boundary condition which applies to this first 
degree differential equation depends on the direction the radiation is traveling. Thus we denote the 
solutions to (103) as /,^(r, +/i) and I„{t, — /i) for upwelhng and downwelhng radiances, respec- 
tively. We shall assume that the upwelling intensity at the surface, Iv{t*, +yu), and the downwelling 
intensity at the top of the atmosphere, 1^(0, — /i), are known quantities. Since fi is positive-definite 
(54), +yU and — /i uniquely specify the angles for which these boundary conditions apply. 

The solution for upwelling radiance is obtained by replacing u in (103) by — /x because u < 
for upwelling intensities. 

^ , = -fi-'e-^/^S+ (104) 

dr 

We could have arrived at (104) by starting from (100b), and proceeding as above except using 
g-T/M as the integrating factor. We now integrate from the surface to level r (i.e., along a path of 
decreasing r') and apply the boundary condition at the surface 

-"/^/,(r,+^) = e-^*/'^/,(r*,+/i)+^-i^ e-^'/^S+ dr' 
/,(r,+/i) = e("-^*)/'^/,(r*,+/i)+/i-i / e("-"')/'^5+ dr' 



e 
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Note that r* > r and r' > r so that both of the transmission factors reduce a beam's intensity 
between its source (at r* or r') and where it is measured (at r). The physical meaning of the 
transmission factors is more clear if we write all transmission factors as negative exponentials. 



(105) 



The first term on the RHS is the contribution of the boundary (e.g., Earth's surface) to the upwelling 
intensity at level r. This contribution is attenuated by the optical path of the radiation between the 
ground and level r. The second term on the RHS is the contribution of the atmosphere to the 
upwelling intensity at level r. The net upward emission of each parcel of air between the surface 
and level r is Sj^{t'), but this internally emitted radiation is attenuated along the slant path between 
r' and r. The factor in front of the integral accounts for the slant path of the emitting mass in 
the atmosphere. 

The solution for downwelling radiance is obtained by replacing u in (103) by /i because fi = 
u in the downwelling hemisphere. The resulting expression must be integrated from the upper 
boundary down to level r, and a boundary condition applied at the top. 





d(e-/^/-) 


Jt'=0 






dr 
t', -fi) 


t'=t 
t'=0 






- 4(0, 


-^^) 


Jo 




e 


^/^4(r, 


-fi) 


= 4(0,-^) + ^-^ Te^'/M^; 

Jo 


dr' 




Ur, 


-fi) 


= 6-^/^(0, -/i)+/i-i / e( 

Jo 










= e-"/^4(0,-/i) / e- 





The upwelling and downwelling intensities in a stratified atmosphere are fully described by 
(105) and (106). Such formal solutions to the equation of radiative transfer are of great heuristic 
value but limited practical use until the source function is known. Note that we have assumed a 
source function and boundary conditions which are azimuthally independent, but that the derivation 
of (105) and (106) does not rely on this assumption. It is straightforward to relax this assumption 
and replace 4(r, /i) by 4(r, /i, 0) and S,y{T, fi) by Si,{t, fi, 0) in the above. 



2.2.2 Thermal Radiation In A Stratified Atmosphere 

Consider a purely absorbing, stratified atmosphere in thermodynamic equilibrium. Then the source 
function equals the Planck function S^j = By = B„{T) (41) and the radiative transfer equation is 
given by (101). It is important to remember that Bi, is the complete source function only be- 
cause we are explicitly neglecting all scattering processes. Thus we need only define the boundary 
conditions in order to use (105) and (106) to fully specify 4. We assume that the surface emits 
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blackbody radiation into the upper hemisphere 

I,{T*,+fi) = B4T{T*)] (107) 

For brevity we shall define 5* = B,y[T{T*)]. At the top of the atmosphere, we assume there is no 
downwelling thermal radiation. 

/,(0, -fi) = (108) 

The solutions for upwelling and downwelling intensities are then obtained by using Su{t', fi) = 
B^{t') (the Planck function is isotropic) and substituting (107) and (108) into (105) and (106), 
respectively 

/,(r,+/i) = e-(^*-^)/'^S,(r*)+/i-i^ e-(^'-^)/^fi,(r') dr' (109) 
4(r,-/i) = /i-i re-(--^')/^i?,(rOdr' (110) 



The first term on the RHS of (109) is the thermal radiation emitted by the surface, attenuated by 
absorption in the atmosphere until it contributes to the upwelling intensity at level r. The second 
term on the RHS contains the upwelling intensity arriving at r contributed from the attenuated 
atmospheric thermal emission from each parcel between the surface and r. The factor in 
front of the integral accounts for the slant path of the thermally emitting atmosphere. The RHS of 
(110) is similar but contains no boundary contribution since the vacuum above the atmosphere is 
assumed to emit no thermal radiation. The upwelling and downwelling intensities in a stratified, 
thermal atmosphere are fully described by (109) and (110). 

2.2.3 Angular Integration 

Once the solutions for the hemispheric intensities are known, it is straightforward to obtain the 
hemispheric fluxes by performing the angular integrations (56a)-(56b). 

F^{T) = 2nr /,(r,-/i)/id/i (111) 

Consider first the downwelling flux in a non- scattering, thermal, isotropic, stratified atmosphere 
obtained by substituting (110) into (111) and interchanging the order of integration 



f-{t) = 2% r ifi ' r \-^^-^'^/^B,{T')dT']fidfi 

J ^=0 y Jr'=0 J 

= 2n / e-("-"')/^5,(r')d/idr' 

Jt'=0 Ju=0 



t'=0 J 11=0 

= 2% [ B^{t')( r e-("-"')/^d/i^ dr' (112) 

Notice that two factors of /i cancelled each other out: The reduction in irradiance due to non- 
normal incidence (fi) exactly compensates the increased irradiance due to emission by the entire 
slant column which is fi^^ times greater than emission from a vertical column. 
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In terms of exponential integrals defined in Appendix 10.7, the inner integral in parentheses in 
(112) isE2(r-r') (603c). 



F~(t) = 2n 



BjT')E,(T-T')dT' 



(113) 



t'=0 



Similar terms arise when we consider the horizontal upwelling flux obtained by substituting (109) 
into (56a) and we obtain 

F+(r) = 271 / /,(r,+^)/id/i 

JtJ.=0 

= 2n (^e-("*-")/'^S^(r*) + /i-i ^ e-("'-")/'^B^(r')dr'^ /id/i 

e-("*-")/>d^ + 27r /" B^{T')(r e'^^'-^^/'^d/i) dr' 

'/i=0 Jt \Jfi=0 J 

= 2nB,{T*)EsiT* -T) + 2n 5,(r')^2(r' - r) dr' (114) 

Jt'=t* 

Subtracting (113) from (114) we obtain the net flux at any layer in a non-scattering, thermal, 
stratified atmosphere 



F+(t)-F-(t) 



2n 



B,{T*)E;ir* -r) + J^^ B,ir')E,ir' - r) dr' - i?,(r')i?2(r - r') d(7l'j5) 



Equations (113) and (114) may not seem useful at this point but their utility becomes apparent in 
§2.3.4 where we define the flux transmissivity. 



2.2.4 Thermal Irradiance 

Assume a non-scattering planetary surface at temperature T emits blackbody radiation such that 
Iv{t* ■, +lA = By(T) (107). What is the total upwelling thermal irradiance from the surface? From 
(56a) we have 



271 ! B^{T)fidn 
Jo 



7t 



2%B,{T) 

2715, (T) 

nB,{T) 

BAT) 



pi 






Jo 






1 


2 





a- 


-0^ 



(116) 
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Notice the isotropy of the Planck function allows the factor of 2 from the azimuthal integration to 
cancel the mean value of the cosine weighting function over a hemisphere. Moving the remaining 
factor of n from the azimuthal integration to the LHS conveniently sets the RHS equal to the Planck 
function. 

We integrate (116) over frequency to obtain the total up welling thermal irradiance 



Thus the factor of it from the azimuthal integration nicely cancels the factor of n from the Stefan- 
Boltzmann Law. Equation (117) applies to any surface whose emissivity is 1. Consider, e.g., a 
thick cloud with cloud base and cloud top temperatures ^(^b) = and T(z^) = Tx, respectively. 
Then the upwelling thermal flux at cloud top and the downwelling flux at cloud bottom will be 

F-{zb) = aT^ and F+(zt) = (tT^, respectively. 

2.2.5 Grey Atmosphere 

Consider an atmosphere transparent to solar radiation and partially opaque to thermal radiation 
governed by (101). The exact solution, including angular dependence, is given in §2.2.2. The hemi- 
spheric fluxes and net flux may only be obtained exactly by accounting for the angular dependence 
as in §2.2.3. We may eliminate the angular dependence of the net flux by making the simplifying 
assumption that the hemispheric up and downwelling irradiances equal a constant times the cor- 
responding intensity. A number of methods exist to determine this constant, called the diffusivity 
factor (e.g., §2.2.15). These methods are all related to the two-stream approximation. 

One such method [Salby, 1996, p. 232] is to identify an effective inclination /2 along which all 
radiation is assumed to travel. With this assumption, the contribution to upwelling irradiance from 
the lower boundary, the first term on the RHS of 115), is 



Inspection (or differentiation) shows that the atmosphere within one optical depth makes the dom- 
inant contribution to (1 18). For Ar = r* — r = 1, the diffusivity factor 




(117) 




(118) 



/i 



-1 



5 

3 



(119) 



The hypothetical coUimated beam of radiation is inclined to the zenith by arccos(3/5) ~ 53.13°. 
This is equivalent to a coUimated beam of radiation travelling vertically through an optical depth 
equal to five thirds the vertical optical depth traversed by the diffuse radiation. 
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Using this assumption (118), the upwelling hemispheric irradiance (56a) for blackbody radia- 
tion is 



2% / I^{+fi)fidfi 
Jo 

2k / /,,(+/i)/id/i 
Jo 

2Kfj,I^{+Ji) / /id/i 
Jo 

',,2" 

2Tlh{+fi) 



2nh{+Ji) 





-0 



(120) 



An analogous relationship holds for the downwelling irradiance. 

Based on (119) and (120), the irradiance structure of the thermal atmosphere may approximated 
by performing a direct angular integration of (100). With our approximation, radiances / integrate 
directly to irradiances F modulo the diffusivity factor ft. (101) 



-fx 



dF- 

dr 
dF+ 



F~ - TiF" 

V V 



/i- 



dr 



B 



(121a) 
(121b) 



It is instructive to examine an idealized grey atmosphere, where the fluxes of interest have 
no spectral dependence. Although this is far from true in Earth's atmosphere, the solution is 
straightforward and sheds light on the radiative equilibrium temperature profile and the greenhouse 
effect. We simplify (121) in two ways. First, we introduce a scaled optical depth df 



I dr. Second, we drop the frequency dependence, which is equivalent to integrating over a broad 
range of frequencies. For heuristic purposes, think of this integration as being over the relatively 
narrow A = 5-20 |im range where most of Earth's terrestrial radiative energy resides. 



dF- 
dF+ 



F- 
F' 



-TlB 
-TlB 



(122a) 
(122b) 



In the absence of dynamical, chemical, and latent heating, the energy deposition in a parcel of 
air is entirely radiative. Under these conditions the idealize grey atmosphere described by (122) 
will adjust to a temperature profile determined by radiative equilibrium. Let the time rate of change 
of temperature T of a parcel be denoted by /i [K s^^], the parcel warming rate^ or cooling rate. The 



^Conventionally h is called the heating rate, which is a misnomer since the dimensions of dT/dt are Ks^^. Thus 
we use the less common but more accurate terminology "warming rate". We reserve "heating rate" for a measure 
power dissipation, i.e., energy per unit time, in J s^^, J m^'^ s^^, or J kg^^ s^^. 
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warming rate is the rate of net energy deposition divided by the specific heat capacity at constant 
pressure Cp [Jkg~^ K~^] times the density p [kgm~^] 

/.^^ = -^^ (123) 

-1 T 



dT 


1 dF 




pCp dz 


K 

s 


V m'^ 



1 



kg K y m s m 

The forcing term on the RHS, dF/dz is the radiative flux divergence, the vertical gradient of net 
radiative flux. Absorption and emission are the only mechanisms which contribute to the flux 
divergence. In terms of hemispheric fluxes, 

^ = ^iF--F^) (124) 
dz dz 

By definition, the time variation of net radiative heating vanishes (dT/dt = 0) at all levels of 
an atmosphere in radiative equilibrium. By (124), we see that the vertical gradient in net radiative 
flux must also vanish in radiative equilibrium. Setting (124) to zero and integrating we obtain 



(125) 



The net radiative flux F(z) = Fq is constant in radiative equilibrium. 
Adding and subtracting (122), we obtain 



4ziF+ - F-) = F+ + F- -2nB (126a) 
dr 

^(F+ + F-) = F+ - F- = Fo (126b) 
dr 

Defining z/^ = F+ - F" and = F = F+ - F", 

— = ^ - 2nB (127a) 
dr 

dib 

= = Fo (127b) 

dr 

Since (j) = — F~ is constant, (127b) shows that d(f)/df = and thus 

ip = 2%B (128) 



Substituting this into (127b), 



^^{2nB) = Fo 

dB _ Fo 

'df ~ 2% 

B{t) = ^ + C, (129) 
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We evaluate the constant of integration by using the boundary condition at the top of the atmos- 
phere. By definition = and f = at TOA. This implies that (p = 'ip{0) = F+(0) = Fq. We 
must therefore have -B(O) = (p/2n by (128). Using this result in (133), 



Substituting (130) back into (133) shows 



271 2% 

271 271 



(130) 



^ ' 271 271 

= §(r + l) (131) 

The thermal absorption and emission in a grey atmosphere increase linearly with optical depth 
from TOA to the surface. 

The upwelling irradiance at the surface is F+(f*) = nB{Ts) where Tg [K] is the surface skin 
temperature. The atmospheric temperature just above the surface is given by (131) as B(f*) = 
|^(f* + 1). Thus there is a temperature discontinuity between the near-surface air and the ground. 

Climate models typically express h (124) in terms of the flux gradient with respect to pressure 
by invoking the hydrostatic equilibrium condition (456) 

^ ^ dT _ 1 dF 
~ dt pcp [-{pg)-^dp] 

Cp dp 

9 (K - F'^) - (-^.+1 - 

Cp Pk - Pk+l 

9 (Fk - Fk^i) + (Fk^i - F^) 
Cp Pk+i - Pk 

where the subscript denotes the A;th vertical interface level in the atmosphere. 



(132) 



(133) 



2.2.6 Scattering 

Energy interacting with matter undergoes one of two processes, scattering or absorption. Scatter- 
ing occurs when a photon reflects off matter without absorption. The direction of the photon after 
the interaction is usually not the same as the incoming direction. The case where the scattered 
photons are homogeneously distributed throughout all An steradians is called isotropic scattering. 
In general, the angular dependence of the scattering is described by the phase function of the inter- 
action. The phase function p is closely related to the probability that photons incoming from the 
direction Q' = {9', (p') will (if scattered) scatter into outgoing direction Cl = {9, (j)). It is usually 
assumed that p depends only on the scattering angle 9 between incident and emergent directions. 

cose = n' n (i34) 
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The case where incident and emergent directions are equal, i.e., Q' = Q corresponds to 6 = 0. 
When the scattered direction continues moving in the forward hemisphere (relative to the plane 
defined by Q'), it is called forward scattering, and corresponds to 6 < n/2. When scattered 
radiation has been reflected back into the hemisphere from whence it arrived, it is called back 
scattering, and corresponds to > n/2. The case where incident and emergent directions are 
opposite, i.e., Q' = —Cl, corresponds to 6 = 71. 

The Cartesian components of ^l' and are straigtforward to obtain in spherical polar coordi- 
nates. 

Cl = sin0cos0i + sin6'sin0j + cos6'k (135) 
O' = sin6''cos0'i + sin6''sin0'j + cos6''k (136) 

The scattering angle 6 is simply related to the inner product of Cl' and Cl by the cosine law 

cose = Cl'Cl 

= sin 6' cos (p' sin 6 cos (f) + sin 6' sin (j)' sin 6 sin (j) + cos 6' + cos 6 

= sin 6' sin 6 (cos 0' cos + sin 0' sin 0) + cos 6' cos 6 

= sm9' sm9cos{(f)' - (p) + cos9' cos9 (137) 



2.2.7 Phase Function 

Accurate treatment of the angular scattering of radiation, i.e., the phase function, is, perhaps, makes 
rigorous demands of radiative transfer applications. A correspondingly large body of literature is 
devoted to this topic. Essential references include van de Hulst [1957], Joseph et al. [1976], M^- 
combe and Grams [1976], Wiscombe [1977], Wiscombe [1979, edited/revised 1996], and Boucher 
[1998]. 

The phase function p(cos 6) is normalized so that the total probability of scattering is unity 

[ p(cose)dfi = 1 (138a) 

1 i>(j>=2Ti i>6=n 

— / j9(^',0';0,0) sin^d^d0 = 1 (138b) 

471 J^=o Je=o 

The dimensions of the phase function are somewhat ambiguous. If the (47r)^^ factor in (138) is 
assumed to be steradians, then p is a true probability and is dimensionless. However, if the (47l)~^ 
factor is considered to be numeric and dimensionless (i.e., a probability), then p has units of (di- 
mensionless) probability per (dimensional) steradian, sr^^. The latter convention best expresses 
the physical meaning of the phase function and is adopted in this text. It is therefore important 
to remember that factors of (47r)^^ which multiply the scattering integral in the radiative transfer 
equation are considered to be dimensionless in the formulations which follow, e.g., (152). Further- 
more, the units of p are probability per steradian, sr~^. 

Scattering may depend on the absolute directions Cl' and Cl themselves, rather than just their 
relative orientations as measured by the angle 6 between them. This might be the case, for exam- 
ple, in a broken sea-ice field. For the time being, however, we shall assume that the phase function 
depends only on 9. 
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In atmospheric problems, the phase function may often be independent of the azimuthal angle 
(j), and depend only on In this case the phase function normalization (138b) simplifies to 

- / p{9';9) sm9d9 = 1 (139) 

2 J 9=0 

In accord with the discussion above, the factor of 1/2 is dimensionless, as is the RHS. 
2.2.8 Legendre Basis Functions 

The phase function (138) specifies the angular distribution of scattering. Solution of any radiative 
transfer equation involving scattering depends on it. Numerical approaches aim for a suitable 
approximation or discretization which represents p(6) to some desired level of accuracy. The 
optimal basis functions for representing p(0) are the Legendre polynomials. 
The Legendre polynomial expansion of the phase function is 

n=N 

p(cose) = ^(2n + l)x„Pn(cose) (140) 

n=0 

An expansion of order contains + 1 terms. The zeroth order polynomial and coefficient are 
identically 1. 

The Legendre polynomials are orthonormal on the interval [—1,1]. The factor of 2r2 + 1 that 
appears in the numerator of the Legendre expansion (146) also appears in the denominator of the 
Legendre polynomial orthonormality property: 

1 /■"=^ 1 

-/ Pn{u)Pk{u)du = ——6n,k (141a) 
^ Ju=-i zn + i 



2 

p„,(cose)Pfc(cose)d(cose) = -— (uib) 

^ Ju=-i Zn + i 

Some other properties of Legendre polynomials are discussed in §10.5. 

The expansion coefficients Xn are defined by the projection of the corresponding Legendre 
polynomial onto the phase function. 



Xn = -f Pn{u)p{u)du (142a) 



u=-l 

e=o 



= -/ P„(cose)p(cose)d(cos0) (142b) 
2 Je=jt 

= -/ P„(cose)p(cose)sinede (142c) 

map the phase function into a series of Legendre polynomials. The Xn are also called the phase 
function moments. Radiative transfer programs like DISORT [Stamnes etal, 1988] usually require 
that directional information be specified as a Legendre polynomial expansion. 

The Legendre expansion for simple, smoothly varying and symmetric phase functions is highly 
accurate with just few terms. For instance, a single moment Legendre expansion exactly describes 
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Rayleigh scattering. However, the expansions of asymmetric phase functions converge much more 
slowly. The strongly peaked forward scattering lobes which appear at large size parameter may 
require hundreds of moments for accurate representation. 

Computing the Xn using quadrature methods significantly increases accuracy and reduces time. 
Wiscombe [1977] recommends Gauss-Lobatto quadrature for highly asymmetric phase functions: 

^ l=L 

Xn = -J2HiPn{cosQi)p{cosQi) sin Qi (143) 
1=1 

where L is the number of Gauss-Lobatto abscissae. Lobatto quadrature is discussed more thorought 
in § 10.6. 



2.2.9 Henyey-Greenstein Approximation 

Computation of the exact phase function of scatterers is laborious but desirable when the objective 
is to predict directional radiances. However, phase function approximations often suffice when 
hemispheric fluxes are the objective. The Henyey-Greenstein Phase Function approximates the 
full phase function in terms of its first Legendre moment, i.e., its asymmetry parameter g. 

= 2^V' nW2 (144) 

[1 + g'^ — 2g cos B)'^/"^ 

The n'th moment in the Legendre expansion of (144) is, conveniently, ^f". 

Xn = 9^ (145) 

Applying (145) in (146) gives 

n=N 

p(cose) = ^(2n+ l)^"P„(cose) (146) 

n=0 

This convenient property (145) makes numerical quadrature of the Legendre expansion coefficients 
unnecessary once the first moment, (? = xi, is known. 



2.2.10 Direct and Diffuse Components 

When working with half -range intensities it is convenient to decompose the downwelling intensity 
I into the sum of a direct component, I~{t, tt), and a diffuse component, I'^{t, tt) such that 

4-(r,i7) = /-(r,i7) + /d(r,r2) (147) 

where we have suppressed the u subscript on the RHS to simplify notation. The direct component 
refers to any photons contributing from a coUimated source which have not (yet) been scattered. 
Typically the coUimated source is solar radiation, and so we subscript the direct component with s 
for "solar". There may be a corresponding direct component of intensity in the upwelling direction 
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/+ if the reflectance at the lower surface is specular, e.g., ocean glint. In such situations it is 
straightforward to define 

/+(r,i7) = /+(r,i7) + /+(r,i7) (148) 

Of course there exist planetary atmospheres somewhere which are illuminated by multiple stars. 
We shall neglect upwelling solar beams for the time being. For consistency, though, we shall shall 
use rather than /+ in equations in which also appears. 

It is plain that (r, tt) is zero in all directions except that of the coUimated beam. More- 
over, the intensity in the direct beam is, by definition, subject only to extinction by Bougher's law 
(i.e., there is no emission). Thus the direct beam component of the downwelling intensity is the 
irradiance incident at the top of the atmosphere, attenuated by the extinction law (91), 

/-(r,i7) = FQe--/'^«5(i7 - rie) (149) 



B.{T) 



2.2.11 Source Function 

The source function for thermal emission S"* is 

k{u) — a{u) 

= {l-w)B,{T) (150) 

With knowledge of the phase function p{tl\ tl) as well as the scattering and absorption coefficients 
of the particular extinction process, we can determine the contribution of scattering to the total 
source function S^. Consider the change in radiance dl^itt) occuring over a small change in path 
df due to a scattering process with phase function p{tl' , tl). The scattering contribution to 
from every incident direction tt' is proportional the radiance of the incident beam, the 
probability that extinction of at location f is due to scattering (not to absorption), and to the 

normalized phase function , (Vj/An (138). 

"^^""^^^ lid') ^^''^ ^^^''^^ 



dr cr(z/) + a{v) 471 

= ^h{n')p{h',n) (151) 

471 

The source function for scattering SI is obtained by integrating (151) over all possible incident 
directions Cl' that contribute the exiting radiance in direction tt 

s: = ^ [ W,n')p{f,Ci',Ci)dn' (152) 
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The LHS and RHS of (152) must have equal dimensions. This is easily verified: tu is dimen- 
sionless, (47r)~^ is dimensionless (cf. §2.2.7), Ii, and SI both have dimensions of intensity, the 
dimensions of p are sr^^ and the cancel the dimensions of dfi' which are sr. 

All terms in (152) are functions of position and scattering process. If there are n distinct 
scattering processes (e.g., Rayleigh scattering, aerosol scattering, cloud scattering) then we must 
know the properties of each individual process (e.g., ai, pi) and sum their contributions to obtain 
to total scattering source function = Xli=i ^ti- neglect such details for now and merely 
remind the reader that applications need to consider multiple extinction processes at the same time. 

The total source function is obtained by adding the source functions for thermal emission (150) 
and for scattering (152) 

s, = {i-w)B,{T) + ^ [ i,if,Ci')pif,Ci',n)dn' (153) 

Although (153) is complete for sources within a medium, additional terms may need to be added to 
account for boundary sources, such as surface reflection. Boundary sources are discussed in §2.3. 



2.2.12 Radiative Transfer Equation in Slab Geometry 

Inserting (153) into (96) we obtain the radiative transfer equation including absorption, thermal 
emission and scattering 

^ - -I +S 
dr 

= -i, + ii-w)B,{T) + ^ [ uif,Ci')pif,n',n)dn' (i54) 

This is a general form for the equation of radiative transfer in one dimension, and the jumping 
off point for our discussion of various solution techniques. The physics of thermal emission, 
absorption, and scattering are all embodied in (154). The most appropriate solution technique 
for (154) depends on the boundary conditions of the particular problem, the fields required in 
the solution (e.g., if irradiance is required but radiance is not), and the required accuracy of the 
solution. 

We are most interested in solving the radiative transfer equations in a slab geometry. To do this 
we transform from the path optical depth coordinate f to the vertical optical depth coordinate r. 
Applying the procedure described in (96)-(100) to (154) we obtain the radiative transfer equation 
for the half range intensities in a slab geometry 

= r(r,{l) - (1 - w)B - — [ I+(T,{l')p(+{l',-{l)dQ' 
dr 471 Jj^ 

f r{T,Ci')p{-Ci',-Ci)dn' (155a) 

47r y_ 

^ d/+(r, n) ^ w)B - — [ i+(T,n')p(+n',+Ci) dn' 

dr 471 

j r{T,Q,')p{-(l',+Q,)dVt' (155b) 



VJ 

471 
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If we substitute (147) into (155a) we obtain 

d/,-(r,r2) d/-(r,r2) _ 



dr dr 



/, (r, 17) + /-(r, Cl)-{l-uj)B-^j^ I;{t, h')p{-h' , -Cl) dQ' 

The direct beam (150) satisfies the extinction law (91) so that — d/^r/dr = . Thus the second 
terms on the LHS and the RHS of (156) cancel eachother and we are left with 

_^ d/,^r,0) ^ Cl)-{1- uj)B - 5*(r, -fl) 

where S* is called the single- scattering source function. S* is defined by the integral containing 
the direct beam on the RHS of (156) 

J i:iT,n')pi-n',-n) dn' 

F0e-^/X-i7Q,-^2) (157) 



w 

471 

w 

471 



Note that (156) is an equation for the diffuse downwelling radiance, not for the total downwelling 
radiance. The relation between the diffuse and total radiance is given by (147). The direct compo- 
nent is always known (150) once the optical depth has been determined. 

Likewise, inserting (147) and (150) into (155b), we obtain the equation for the diffuse up- 
welling intensity 

^ d/+(7-,17) ^ zu)B - 5*(r, +n) 



dr 

I 

471 



^ / i+{T,n')p{+n',+n)d9! I q{T,n')p{-n',+n)d9! (i58) 

471 J , 471 J_ 



where 



s*{r,+n) = ^Ji;{T,n')p{-n',+n)dn' 

j FQe-^^''°5{n',CiQ)p{-Ci',+Ci) dn' 

-FQe-^/X-f2Q,+f2) (159) 



w 

4% 

w 

471 " 
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2.2.13 Azimuthal Mean Radiation Field 

Often we are not concerned with the azimuthal dependence of the radiation field. In some cases 
this is because the azimuthal dependence is very weak. For example, heavily overcast skies, or 
diffuse reflectance from a uniform surface. In other cases the azimuthal dependence may not be 
weak, but there insufficient information to fully determine the solutions. For example, the full 
surface BRDF or the shapes of clouds are not available. The azimuthal dependence of a radiative 
quantity X(0) is removed by applying the azimuthal mean operator 



= — X(0)d0 (160) 

For future reference we present also the polar angle integration operator which will be used to 
formulate the two stream equations. 



Xe = / X(^)sin^d^ 
Jo 

= [ X(/x)d/i (161) 



Applying (160) to (25) and (138), we obtain the azimuthal mean hemispheric intensity, phase 
function, and single-scattering source function respectively. 



/^(r,/i) = /^(r,/x,0)d0 (162) 

p(±/i';±/i) = — / p(±/i',0';±/i,</>)d</) (163) 

Jo 



S*iT,±fx) = £Foe--/X-^o,±/i) (164) 

No additional symbols are used to indicate azimuthal mean quantities. The presence of only the 
polar angle (9 or fi) on the LHS of (162) indicates that azimuthal mean quantities are involved. 

Applying (160) to the full radiative transfer equations for slab geometry, (156) and (158), we 
obtain the radiative transfer equations for the azimuthal mean intensity in slab geometry 

_^d^^^ = /d (r,/.) - (1 - - £Fee-^/X-/zo, -/x) 

~ f X ^d^(^' -f^) d-"' - f / -/^) d/i' (165a) 

^d^^ = I^T, /x) - (1 - w)B - ^Foe--/X-/io, 



2 



J Idi^^ +/^) d/^' ~ ^ J h +/^) d/i' (165b) 



The azimuthal mean equations (165) are very similar to the full equations (156) and (158). While 
the remainder of the terms in (165) contain only one azimuthally dependent intensity, the scattering 
integrals contain two, p and I"^. This causes the factor of zu/ An in front of the scattering integrals 
has to become a factor of ct7/2. 
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We apply the two stream formalism introduced in §2.4 to the radiative transfer equation for the 
azimuthal mean radiation field (165). Recall that the two stream formalism is obtained by applying 
three successive approximations. The approximation is achieved in three stages. First, we operate 
on both sides of (165a) with the hemispheric averaging operator (161). The procedure for (165a) 
is identical. 

- f^^^"^^^ d/i = (r, fi) d^ - (1 - w) j^' Bdfi- ^F^e-^/^' j9(-/io, -/i) d/i 

£ /d^ (r, -/i) d/i' d/i - I (r, fi')p{-f^', -/i) d/i' d/i (166) 



zu 
~2 



Further approximations are required to simplify equations (173). These approximations allow 
us to decouple products of functions with /i dependence. The first approximation, (167), replaces 
the continuous value ji on the LHS of (165) by a suitable hemispheric mean value Jx. 

' d/^-(r,/x)^ d r\_^ d/d(r) 

1 d/i ^ /i— / /^(r,/x)d/i = /i^^ (167) 







dr dr L dr 



Extracting /i from the integral on the LHS of (167) is an approximation whose validity worsens as 
the correlation between fi and /^(r, fi) increases. 

The first two terms on the RHS of (166) are simple hemispheric integrals of hemispherically- 
varying intensities. We will replace /^(r, /i) by the hemispheric mean intensity /^(r) (205). The 
Planck function is isotropic and so may be pulled outside the hemispheric integral which promptly 
vanishes since d/i = 1. Thus the hemispherically integrated intensities which result explicitly 
are replaced by the hemispheric mean intensity (205) associated with jl. 

The final three terms on the RHS of (166) involve products of the phase function and the 
intensity. We approximate the scattering integrals by applying the hemispheric integral over fi to 
the intensities, and then extracting the intensities from original integrals over /i'. 

/ / P(+/^', -/^)^/(^,/^')d/i' ^ / /^(r,/i')d/i / p(+/i', -/i) d/i' 
Jo J+ Jo J+ 

= /d+(r)£p(+/i',-/^)d/i' (168) 

These terms define the backscattered fraction of the radiation field. Equation (168), for example, 
represents the fraction of upwelling energy backscattered into the downwelling direction. This is 
one of four similar terms that appear in the coupled equations (166). We define the azimuthal mean 
backscattering function h{jj) as 

W = \ P{-f^', /^) d/^' = ^ ^ -/^) d/i' (169) 
The complement of the backscattering function is the forward scattering function f. 
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Pairs of the terms in (169)-(170) are identical due to reciprocity relations. Reciprocity relations 
state that photon paths are reversible. 

The hemispheric mean backscattering coefficient b is the hemispheric integral of (169) 

b = f b{fi)d^ (171) 



7) [ [ p(-Ai',/i)d/i'd/i (172) 
^ Jo Jo 



'0 ^0 

The hemispheric mean forward scattering coefficient is defined analogously. 
The result of these two steps is 

-/i^^ = lair) - (1 - ^17)5 - ^Foe--/X-/^o, -/i) 



Y / ^d"(^)p(+/^''-^)d^' - ^ J hir)pi-f^\-l^)dfx' (173a) 



fi^^ = /+(r) - (1 - v.)B - ^F^e'^/^^p{-^,o, 



to 
~2 



£ /+(r)p( V, +/.) d/i' - I /d V, d/i' (173b) 

These are the hemispheric mean, azimuthal mean, two stream equations for the radiation field in 
slab geometry. It should be clear that these equations (173) are not derived simply by performing 
a hemispheric integral (161) on (165). A strict hemispheric averaging operation would applies to 
the product of the phase function and and the intensity, not to each separately, so that extracting 
from the scattering integrals is an approximation, as is the definition of Ji. 

2.2.14 Anisotropic Scattering 

Armed with techniques introduced in solving the two stream equations for an isotropically scatter- 
ing medium §2.4.1, we now solve analytically the coupled two stream equations for an anisotropic 
medium. First we recast (??) into a simpler notation 

d/+fr) 

/i^^ = /+(r) - - b)I^{T) - wbl;{T) - (1 - w)B, - (174a) 
dr 

-f.^LM. = j-^r) - w{l - b)I;{T) - wbI+{T) - (1 - uj)B, - S;- (174b) 
clr 

For the remainder of the derivation we drop the u subscript and the explicit dependence of I"^ on 
T. Adding and subtracting we obtain 

d(/+ + /-) 



dr 
d(/+ - /- 

dr 



-{a- I3){r - r) (175a) 
-(« + /?)(/+ + r) (175b) 



where we have defined 



« = -[1 -6)]//2 (176a) 
(3 = wb/ji (176b) 



2 RADIATIVE TRANSFER EQUATION 



47 



2.2.15 Diffusivity Approximation 

The angular integration required to convert the intensity field into a hemispheric irradiance is a 
time consuming aspect of numerical models and should be avoided where possible. According to 
(56a)-(56b) 

1 

-F+ = 2 / /,(+/i)/id/i (177a) 

7t ./n 

1 



1 

F; = 2 I /,(-/i)/id/i (177b) 



The factor of n has been moved to the LHS so that, when the source function is thermal, the RHS 
is the integral Planck function (116). These angular integrals may be reduced to a calculation 
of the intensity at certain quadrature points (zenith angles) in each hemisphere with surprising 
accuracy. The location of the optimal quadrature points may be arrived at through both theoretical 
and empirical methods. 

Two point /wZ/ Gaussian quadrature (§10.5) tells us the optimal angles to evaluate (177) at are 
u = ±3-1/^ e = ±54.7356°. 

^F- ^ -|/^(_3-i/2) (178b) 

The relation in Equation (178) is exact when luil^) is a polynomial of degree < 2. For two point 
Gaussian quadrature the quadrature weights happen to be unity, but not so for three point (or higher 
order) Gaussian quadrature. Here is the four point quadrature version of (178) 

Ml = -M4 = 0.339981 

U2 = -Us = 0.861136 

Ai = = 0.652145 

A2 = As = 0.347855 (179) 

— F+ AiUih{ui) + A2U2lu{u2) 

2% 

—F- ~ A-iUsh{u-i) + AiuJ^{ui) 
2% 



What is often used to evaluate (177) instead of Gaussian quadrature is a dijfusivity approxi- 
mation. Instead of choosing optimal quadrature angles /i^, the diffusivity approximation relies on 
choosing the optimal effective absorber path. Thus the dijfusivity factor, D, is defined by 

-F+ ~ I,{u = +D-^) (180a) 
It 

-F- - h{u = -D-') (180b) 

Comparing (180) to (178) and (177) shows that D replaces, simultaneously, the quadrature angle, 
its weight, the factor of n in (177), and a factor of 2 from the azimuthal integration. The angle 
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6d = arccosD^^ may be interpreted as the mean slant path of radiation in an isotropic, non- 
scattering atmosphere. 

The reasons for subsuming so many factors into D are historical. Nevertheless, (180) is much 
more common than (178) in the literature. Heating rate calculations in isotropic, non-scattering 
atmospheres show that using D = 1.66 results in errors < 2% [Goody and Yung, 1989, p. 221]. 

In applying (180), D only affects the optical depth factor of the intensity. Thus computing 
the intensity in (180) at the angle 9 = is formally equivalent to computing the intensity in 
the vertical direction = but through an atmosphere in which the absorber densities have been 
increased by a factor of D. 



2.2.16 Transmittance 

The transmittance between two points measures the transparency of the atmosphere between those 
points. The transmittance T from point Pi to point P2 is the likelihood a photon traveling in 
direction Q at Pi will arrive at P2 without having interacted with the matter in between. 

T(Pi,P2)=exp[-f(Pi,P2)] (181) 

We have left implicit the dependence on frequency. In a stratified atmosphere, we are interested in 
the transmission of light traveling from layer z' to layer z at angle 9 from the vertical, thus 

T(z',2;;/i) = exp[-r(z',z)//i] (182) 

Thus, in common with absorptance A and reflectance R, the transmittance T G [0, 1]^ 

The vertical gradient of T is frequently used to formulate solutions of the radiative transfer 
equation in non- scattering, thermal atmospheres. 

d d 

—T{z',z;^m) = —exp[-T{z',z;fi)/fi] 

dT{z',z;fi) 1 r . / w .dT{z',z;iJ,) 

dz' = --^-pM^'^;/^)//^]^^ 

1 dr 

= —T{z',z;fi)— (183) 
fl oz' 

Inverting the above, we obtain 

-irri / N dT{z',z;fi) dz' 
fx T{z,z;ix) = — — (184) 

The integral solutions for the upwelling and downwelling intensities in a non-scattering, ther- 
mal, stratified atmosphere (109)-(1 10) may be rewritten in terms of T. (110) from r to z. To do 



^This text distinguishes between transmittance and transmission. The former measures the potential for the latter to 
occur. Thus a transmittance of unity means the transmission of energy is very high. A similar rule holds for reflectance 
and reflection and for absorptance and absorption. 
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this we use 



z{t = T*) 


= 


z{t = 0) 


= oo 


T(0,2;;/i) 


— g-(T*--r)/M 


T{z',z;n) 




T{z,z';n) 




dT{z',z;fi) 


dT{z, z' 


dz' 


dz' 



(185) 

The last relation simply states that the transmission decreases as the distance between z! ^ z in- 
creases and visa versa. Substituting the above into (109) and (110) yields 

Iu{z,+fi) = T{0,z;fi)B,{0)+fi-^ J T{z',z;fi)B,{z'){-dz') 
h{z,-ii) = fi-' [ T{z,z';fi)B,{z')i-dz') 

J oo 

Rearranging terms we obtain 

= T(0,2;/i)5,(0)+/i-^ /" T{z',z;^)B,{z')dz' (186a) 

/•oo 

I,{z, -fi) = fi-' / T{z, z'- fi)B,(z') dz' (186b) 

J z 

Substituting (184) in the above leads to 

I,{z, +fi) = T(0, z; fi)B,{0) + ^I^^^B^iz') dz' (187a) 

I,{z, -/x) = - ^lih^lAB^iz') dz' (187b) 

If T is known, then Equations (187a)-(187b) are suitable for use in quadrature formulae to obtain 
irradiances. 



2.3 Reflection, Transmission, Absorption 

Perhaps the most useful metrics of the radiative properties of entire systems are the quantities re- 
flection, transmission, and absorption. These metrics describe normalized properties of a system 
and thus are somewhat more fundamental than absolute quantities (like transmitted irradiance), 
which may, for example, change depending on time of day. We shall define the reflection, trans- 
mission, and absorption first of the radiance field, and then integrate these definitions with the 
appropriate weighting to obtain the R, T, and A pertinent to fluxes. We shall introduce a painful 
but necessary menagerie of terminology to describe the various species of R, T, and A. 

Molotch et al. [2004] show how remotely sensed surface reflectance differs significantly from 
that derived from snow-age models. Using the more realistic albedos improved models of snowmelt 
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by providing more accurate net surface shortwave radiation estimates. This term, dominates 
the snow-melt energy budget. 

Wang et al. [2005] derived a simple functional fit to the desert surface albedo observed by 
MODIS. 

2.3.1 BRDF 

The bidirectional reflectance distribution function (BRDF) p is the ratio of the reflected intensity to 
the energy in the incident beam. As such, p is a function of frequency, incident angle, and scattered 
angle 

piv,-h',h) = — — (188) 
i-{n')cose'dQ' 

The dimensions of p are sr~^ and they convert irradiance to intensity. The reflected intensity 
in any particular direction fl is the sum of contributions from all incident directions Cl' that have a 
finite probability of reflecting into Cl. 



= f ^i;{Ci')p{u,-Ci',Ci) cose' an' (i89) 

The dependence of p(z/, —Cl', Cl) on two directions and on frequency makes it a difficult property 
to measure. 

2.3.2 Lambertian Surfaces 

Fortunately many surfaces found in nature obey simpler reflectance properties first characterized 
by Lambert. A Lambertian surface is one whose reflectance is independent of both incident and 
reflected directions. The reflectance of a Lambertian surface depends only on frequency 

p{u,-n',n) = pj^iu) (190) 

The intensity reflected from a Lambertian surface is given by inserting (190) into (188) and then 
using (56b) 



i^rW = f i;mpM cose' dn' 
= pl(z^) f .^i; (Ci') cose' an' 



= pMF- (191) 

As expected, the intensity reflected from a Lambertian surface depends only on the incident irradi- 
ance, and not at all on the details of the angular distribution of the incident intensity field. 
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The irradiance reflected from a Lambertian surface Fj^ is the cosine-weighted integral of the 
reflected intensity (191) over all reflected angles 



FX = L /+(^i)cos0d^] 



Pi^{u)F^ cos edVt 
1~2 



Neither F^ and Pl(^) depend on the emergent angle VI so the integral reduces to the familiar 
integral of cos over the hemisphere (59) which is n. 

F+ = npi.{u)F; (192) 

The reflected irradiance may not exceed the incident irradiance or the requirement of energy 
conservation will be violated. Therefore (192) shows that 

< T^-' (193) 

with equality holding only for a perfectly reflective (non-absorbing) Lambertian surface. 

Many researchers prefer the Lambertian BRDF to have an upper limit of 1, not 71^ (193). Thus 
it is common to encounter in the literature BRDF defined as r = np. 

2.3.3 Albedo 

The reflectance of a surface illuminated by a coUimated source such as the Sun or a laser is of great 
interest in planetary studies and in remote sensing. For an incident coUimated beam of intensity 
Fq, the diffusely reflected intensity is 

/+ = / ^ Fq5{(1' - CIq)p{u, -Cl', Cl) cos 9' dfi' 

= Fqp{u,-{Iq,{1)cos9o (194) 
Integrating (194) over all reflected angles we obtain the diffusely reflected hemispheric irradiance 

F+ = I FQp{iy,-nQ, n) cos 9o cos 9 dn 

J+Cl 

= Fqcos9o / . p{u,-CIq,Ci) cos 9 dQ (195) 



The flux reflectance or plane albedo is the ratio of the reflected (diffuse) irradiance (195) to the 
incident (coUimated) irradiance (150) 

p(z.,-f2Q,27r) = ^ 

= I ^ p{u,-{lQ,{l)cos9dQ (196) 

The 271 indicates that the plane albedo pertains to the entire reflected hemisphere. 
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Table 4: Surface Albedo" 

Surface Type Reflectance 



Glacier 


0, 


.9 




Snow 


0, 


.8 




Sea-ice 


0. 


.D 




Clouds 


0, 


.2- 


0.7 


Desert 


0, 


.2- 


0.3 


Savannah 


0, 


.2- 


0.25 


Forest 


0, 


.05 


-0.1 


Ocean 


0, 


.05 


-0.1 


Planetary 


0.3 





"Sources: Bird Entrails (2000) 

When the reflecting body is of finite size then the corresponding ratio of reflected to incident 
fluxes is more complicated because it contains "edge effects", e.g., diminishing contributions from 
planetary limbs. First we consider the reflection, called the spherical albedo, planetary albedo or 
Bond albedo, of an entire planetary disk illuminated by coUimated sunlight. 

The directional reflectance p(z/, —271, Cl) is the 

= f^F-{-{l')p{u,-2Ti,n)cos9'dQ' 

(197) 

The flux reflectance of a Lambertian surface to coUimated light is found by subsituting (190) 
in (196) 

— i^Q, 271) = / ^ pl('^) cos 6* df2 

J +d 

= TipM (198) 

which agrees with (192). 

The BRDF must be integrated over all possible angles of reflectance in order to obtain the flux 
reflectance, p{u, — f^', 27r). Representative flux reflectances of various surfaces in the Earth system 
are presented in Table 4. 



2.3.4 Flux Transmission 

T\\Qflux transmission between two layers is (twice) the cosine- weighted integral of the spectral 
transmission 

Jo 
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T,^{z',z) = 2 [ exp[-r{z',z)/fi]fidfi (199) 

We can rewrite (199) in terms of exponential integrals (603d) 

T^^{z',z) = 2E3[t{z',z)] (200) 

The factor of 2 ensures that expressions for the hemispheric fluxes in non- scattering, thermal at- 
mospheres closely resemble the equations for intensities, as will be shown next. 

The vertical gradient of the flux transmission is frequently used to formulate solutions of the 
radiative transfer equation in non- scattering, thermal atmospheres. It is obtained by differentiating 
(200) then applying (602a) 

-T^{z,z) = 2—Es[t{z,z)]- 



dz' ^ ' dz' ' dz' 

dz' - 2h^2(^,^)J^^ 

dr 

= -2E2{z',z)— (201) 

The integral solutions for the upwelling and downwelling fluxes in a non-scattering, thermal, 
stratified atmosphere (113) and (114) may be rewritten in terms of T^. The procedure for doing 
so is analogous to the procedure applied to intensities in §2.2.16. First, we change the independent 
variable in (1 13) and (1 14) from r to z. 

rz'=z 

F+{z) = 2nB,{0)E;{z) + 2n B,{z')E2{z - z') dz' 

Jz'=0 

/>z'=oo 

F-{z) = 2n B,{z')E2{z' -z)dz' 

J z'=z 

Substituting (201) into the above yields 

F^iz) = 5.(0)Tj^(0, z) + / B^iz' f "^Z ' dz' (202a) 



n ^ A,=o ' dz^ 

1 ^^'=00 dT^iz' z) 

-F; {z) = - B,{z') : ' dz' (202b) 

% Jz'=z 

The analogy between (202a)-(202b) and (187a)-(187b) is complete. plays the same role in 
former that T plays in the latter. 

Equations (202a)-(202b) are particularly useful because they contain no explicit reference to 
the angular integration. The polar integration, however, is still implicit in the definition of T^f 
(199). The hemispherical irradiances may be determined without any angular integration if the 
diffusivity approximation is applied to (202a)-(202b). This is accomplished by replacing by 
the vertical spectral transmission T through a factor D (180) times as much mass 

Tj'(z', z) ~ T(z', z- arccos D'^) (203) 



Determination of T requires the use of a spectral gaseous extinction database in conjunction with 
either a line-by-line model, a narrow band model, or a broadband emissivity approach. 
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2.4 Two-Stream Approximation 

The two-stream approximation to the radiative transfer equation assumes that up- and down- 
welling radiances travel at mean inclinations /2+ and p,^ , respectively. With this assumption, the 
hemispheric intensities if in the azimuthally averaged radiative transfer equation (165) lose their 
explicit dependence on fi and depend solely on optical depth. For an isotropically scattering at- 
mosphere in a slab geometry, the two-stream approximation may be written 

/^^^^ = - |/+(r) - |/;(r) - (1 - w)B^[T{t)] (204a) 



d/+(r) 


= 


dr 




d/-(r) 




dr 





-/i-^^ = lAr) - ^Kir) - ^I-{r) - (1 - w)BAT{r)] (204b) 

These equations are no longer exact, but only approximations. Nevertheless, and as we shall show, 
the two-stream approximation yields remarkably accurate solutions for a wide variety of scenarios 
including the most important ones in planetary atmospheres. We emphasize that the solutions to 
the two-stream equation. If, are hemispheric mean intensities each associated with a hemispheric 
mean direction Ji^. In this interpretation of Ifir), 



Jo d/x ^0 



(205) 



As emphasized by (205), If are not isotropic intensities in each hemisphere. 

We now re-define the other important properties of the radiation field to be consistent with the 
two-stream approximation. Perhaps most important is the hemispheric irradiance. Starting from 
(56) we obtain 



F±(r) = 271 / J±(r,/i)/xd// 

^0 





« 2nlilf{T) (206) 

where ~ indicates that (206) is an approximation. Equation (206) reveals one interesting physical 
interpretation of the two-stream approximation: the radiance field comprises only two discrete 
streams travelling at angles /2^. Mathematically, this radiance field could be represented by the 
union of the two delta-functions representing coUimated beams travelling in the directions. 
This interpretation (presumably) gives the "two stream" method its name. Physically, it is more 
appropriate to interpret the two-stream radiance field as being continuously distributed in each 
hemisphere (205), such that the angular moments (e.g., in the definition of irradiance) have the 
properties in (206). 

A common mis-understanding of the two-stream approximation is that the radiance field is 
isotropic in each hemisphere. If two-stream radiance were isotropic in each hemisphere then the 
appropriate measure of hemispheric irradiance would be Ff = nif (59) as opposed to Ff = 
2n]ji^lf (206). In fact, assuming the intensity field varies linearly with /i (and is thus anisotropic) 
is a reasonable interpretation of the two-stream approximation. 

The mean intensity I^j (20) for an azimuthally independent, isotropically scattering radiation 
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field is 




Iu{t, u) du 



(207) 



Although luij) (207) appears independent of /i, the two-stream intensities in (207) do depend on 
/i. 

Continuing to write the two- stream approximations for radiative quantities of interest, we turn 
now to the source function. The source function Sy (153) for an azimuthally independent, isotrop- 
icaUy scattering (p = 1) radiation field is 



The final term on the RHS is the two-stream mean intensity ntnmnfrq (207). ly appears in 5*,^ 
(208) due to the assumption of isotropic scattering. For more general phase functions, Sy will not 
contain 1^. 

A comment on vertical homogeneity is appropriate here. The Planck function in (208) varies 
continuously with r, whereas w is assumed to be constant in a given layer over which (208) is 
discretized. Thus r-dependence appears explicitly for B^, and is absent for w. More generally, 
the radiation field (intensity, irradiance, source function, etc.) varies continuously with r. The 
discretization of r (or, in general, all spatial coordinates) is required in order to solve for the con- 
tinuously varying radiation field. The discretization determines the scale over which the properties 
of the matter in the medium may change. 

The radiative heating rate in the two-stream approximation may be obtained in at least two 
distinct forms. 



The first term on the RHS represents the rate of absorption of radiative energy, while the second 
term is the rate of radiative emission. Note that (209) depends on the absorption cross section a 
rather that the single scattering albedo vj. This emphasizes that scattering has no direct effect on 
heating. 

Together the relations for Fy, ly, Sy, and hy in (206)-(209) define a complete and self-consistent 
set of radiative properties under the two-stream approximation. We turn now to obtaining the two- 
stream intensities lyi^r) from which these quantities may be computed. 




(l-t^)5.[r(r)] + -[/+(r) + C(r)] 



(208) 



dFy 



dz 

4'!lIy{T)a - AnaBy[T{T)] 
27ra[/+(r) + l;{r)] - 4TiaBy[T{T)] 



(209) 
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2.4.1 Two-Stream Equations 

We shall simplify both the physics and nomenclature of (204) before attempting an analytical 
solution. First, for algebraic simplicity we set /2+ = /i~ = fx. Also, we neglect the thermal source 
term B^[T(t)]. With these assumptions, (204) becomes 



^iim _ J, (r)^^Ii(r)- |/; (r) (2 10a) 



-/i^ = 4-W-f/^(r)-f4-(r) (210b) 

For the remainder of the derivation we drop the u subscript and the explicit dependence of 
on r. Adding and subtracting the equations in (210) we obtain 

^ j+ ^ J- _ ^ j-^ ^ j-^ (211a) 

dr 

i.^AIl±n ^ J, _ J- (211b) 

dr 

These may be rewritten as first order, coupled equations in Y"^ where 

= /+ ± I (212a) 

/± = l(y+±F-) (212b) 



Using (212a) in (211) 

dr- d - zu) 



-r+ (213a) 



dr fj, 
dY+ 1 

= -y- (213b) 

dr fi 

We differentiate (213) with respect to r to uncouple the equations, 

d^y- _ (1 - w) dy+ _ (1 -U7) 1^^, _ 1 - zu ^_ 

dr^ /i dr Ji Ji Ji^ 

dr^ dr Ji Ji Jj? 

The uncoupled quantities satisfy the same second order ordinary differential equation 

d?Y^ 

dT2 

where 



r^F^ (214) 



^ (215) 

Elementary differential equation theory teaches us that the solutions to (214) are 

F±(r) = Cfe^^ + Cte-^^ (216) 
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The general forms of are obtained by inserting (216) into (212b) 

/+(r) = Cie^" + Cae"^" (217a) 
/- (r) = Cge^" + C^e-^^ (217b) 

Equation (217) suggests that there are four unknown constants of integration C1-C4. This is 
a consequence of solving for the hemispheric intensities (25) separately in two equations (210) 
rather solving the full radiative transfer equation (154) for the full domain intensity /. As dis- 
cussed in §2.1.16, the motivation for solving the half -range equations is that the physical boundary 
condition (at least in planetary science applications) is usually known for half the angular domain 
at both the top and the bottom of the atmosphere. These two independent boundary conditions for 
the first order equations for the hemispheric intensities (210) mean there must be two additional 
relationships among the unknowns C1-C4. Note that the boundary conditions themselves are still 
unspecified. 

To obtain the relationships among C1-C4 we substitute the general solution forms (217) into 
the governing equations (210). It can be shown (after much tedious algebra) that 

C2 = P00C4 

C3 = PooCl 

where 



p» . (218) 

With these additional constraints, (217) becomes 

/+(r) = Cie^" + PooC^e-^^ (219a) 
r (r) = PooCie^^ + C^e-^^ (219b) 

where now only the constants Ci and C4 remain to be determined by the particular boundary 
conditions of the problem. 

The most relevant two-stream problem that is easily tractable is that of an isotropic, down- 
welling hemispheric radiation field of intensity X illuminating a homogeneous slab which overlies 
a black surface (so that upwelling intensity at the lower boundary is zero). 

/+(r = r*) = (220a) 
/-(r = 0)=X (220b) 

Substituting (220) into the solutions for the case of an isotropically scattering medium (219) 
we obtain 

/+(r = r*) = = Cie^"* + p^C^Q-^^' 
/-(r = 0) = X = p^Ci + C4 
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which lead to 



= , = -Y^^ (221a) 

Xe'"'^* X 
C4 = -FT* 2 — ^rv^ ^ 1 2 — -2Vt* (221b) 



It will prove convenient to define the denominators of (221b) as 

V = e^"* - p^e"^"* (222a) 
V* = 1- p^e"^^"* (222b) 

where V* = e^^^*V. The initial forms of Ci and C4 contain both positive and negative exponen- 
tials and are pleasingly symmetrical. These forms have been divided top and bottom by e^^""* to 
obtain the final forms on the RHS for reasons that will be discussed shortly. 

Substituting (221) and (222) into (219) we can finally write the intensities at any optical depth 
r in terms of the known quantities r* and zu 

= P^[gr(.*-r) _ e-r(.*-r)] ^ ^[e-r- _ e-^H— /2)] (223a) 

I-{r) = ^[e^r'-r) _ ple-nr^-r)] ^ |_^^-r. _ pl,-mr'-r/2)^^ (223b) 

Two forms of the solution are presented. The initial forms on the RHS are differences between 
positive and negative exponential terms. In this form, the solutions may appear to depend only on 
the distance from the lower boundary (r* — r) but notice that the V term (222a) depends on the 
absolute layer thickness (r*) as well. However these forms are not recommended for computational 
implementation because the positive exponentials are difficult to handle for large optical depths. 
Moreover, the difference between the exponential terms will quickly lead to a loss of numerical 
precision as r — )• oo. 

The final forms on the RHS of (223) contain only negative exponentials and are numerically 
well-behaved as r — )■ oo. Since < r < r*, it is clear that all exponentials in (223) and (222a) 
are negative exponentials and thus well conditioned for computational applications. There are 
no physical processes which would lead to positive exponential terms in the solutions, so it is 
always worthwhile doublechecking for accuracy or stability any formulae which contain positive 
exponentials. 

As expected, the solutions (223) are proportional to the incident intensity X, which is the only 
source of energy in the problem since we neglected the thermal source term while constructing 
the governing equations (210). The upwelling radiance is proportional to the parameter poo- The 
physical meaning of poo is elucidated by noting that /+ — )• poo^^ as r* oo. Thus poo is the 
maximum reflectance of semi-infinite layer of matter with the same optical properties as the layer 
in question. 



2.4.2 Layer Optical Properties 



The solutions for the intensities (223) allow us to derive the other optical properties of the layer 
using (206)-(209). The total flux reflectance p{v, —271, 27r) is the ratio of reflected irradiance to all 
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0) ^, 27r/i/+(r = 0) /J 



r 



0) 



F-(r = 0) 27r/i/-(r = 0) /-(r 

-1 



0) 



V* 

Poo 

V* 

Poo(l 



1 



1 - e 



2rT 



-2rT* 



X 



X 



Poo ^„rT* 



' — e " ' ) 



-rr* 

-^T*^ 



1 - P^e 



-2rT* 



-.Ft* 



P^e-r-* 



(224) 



where we have used the definition of V* (222b) in the last step. 

The total flux transmittance T{u) is the ratio of transmitted irradiance to incident irradiance. ^° 



T{u) 



T{u, -271, -271) 
I 



[r = T 



-Ft* 



V* 



pie-'- 



V* 



2n-pI~{T = T*) _ I-{t = t* 

27r^/-(r = 0) 

-1 

) 



0) 



1 - pLe-^^^ 



(1 - Ppe-"- 
V* 

(1 - Ppe-"" 
1 - p^e-2r^* 



1-p: 



V 



pIc 



(225) 



The first expressions in (224) and (225) were derived using the V* forms of the solutions (223) 
and contain no positive exponential terms which blow up as r* — > oo. The second expressions, on 
the other, hand, are numerically ill-conditioned as r* — ^ oo. 

Conservation of energy requires that the total flux absorptance A{v) is one minus the sum of 
the transmittance and the reflectance 



A{V) = A{U, -2%) : 

= 1-T{u)-R{iy) 



F~iO) - F-{t*) + F+{t*) - F+(0) /-(O) - /-(r*) + /+(r*) - /+(0) 



4-(o) 



(226) 



2.4.3 Conservative Scattering Limit 

The preceding section evaluated -R(z^) and r(z/) for the entire layer for a non-conservative scat- 
tering (zu 7^ 1) atmosphere. When zu = 1 the two-stream equations take a simpler form with the 

^"A non-reflecting lower boundary ensures there is no ambiguity to the meaning of transmittance. A reflecting 
lower boundary artificially enhances the apparent transmittance by allowing surface-reflected radiation to scatter off 
the bottom of the layer and thus be counted as transmitted radiation. 
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J+Cr = 0) 



Ar* - r) 

2fl + T* 

r*-r 

2jl + T* - 

T* 
2J1 + T* 



X 



4-(t = o) 

X(2/2 + T* -t)" 
2/i + T* 



T(.) . T(.,-2.,-2.).|(^ 

X{2li + T*-T) 1 

2/i + T* X 

_ 2/i + T*-T 

~ 2/i + T* 
_ 2ji 

~ 2fl + T* 
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2.5 Example: Solar Heating by Uniformly Mixed Gases 

In hydrostatic balance, the density of a uniformly mixed species i is 

p,{z) = p,(z = 0)e-^/^ (229) 

(230) 

The optical depth of a purely absorbing species with volume absorption coefficient /ca is (84) 

fxm (231) 

2.6 Exercises for Chapter 2 

1. Consider a vertically homogeneous stratiform liquid water cloud extending from z = 1- 
2 km above the surface, i.e., the cloud is 1 km thick. Assume the sun is directly overhead and 
let F{z) denote the downwelling flux (in W m^^) in the direct solar beam. At the top of the 
cloud F{z = 2) = Fq. In the middle of the cloud it is found that F{z = 1.5) = Fq/2. 

(a) What causes F(z) to decrease from the top to the middle of the cloud? 
Answer: Absorption and scattering by cloud droplets. 



(b) What is F{z = 1), i.e., the downwelling flux exiting the bottom of the cloud? Answer: 

^0/4 

(c) What is the extinction optical depth of the cloud r^l Answer: In 2 

(d) Assuming the cloud droplets are uniform in size, what additional information about 
the cloud is needed to estimate the droplet radius? Answer: The mass or number 
concentration of cloud droplets and their density (i.e., the density of water, 1 gcm^^). 
Also, the extinction efficiency Qe must be known, but it can be assumed to be 2. 
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3 Remote Sensing 
3.1 Rayleigh Limit 

The size parameter x is the ratio of particle circumference to wavelength, times the real part of the 
refractive index of the surrounding medium 

X = krn^ (232) 
= 2TirnJX (233) 

For particles in air, ^ 1 so x ^ 2Tir/\. The Rayleigh limit is the regime where the wavelength 
is large compared to the particle size, A ^ r, or, equivalently, where x — )■ 0. Note that this 
regime applies equally to visible light interacting with nanoparticles, and to microwave radiation 
interacting with cloud droplets. The absorption and scattering efficiencies reduce to simple, closed 
form expressions in the Rayleigh limit 

lim^Qa = -4x^>(r^^) (234a) 

(234b) 



''X 



^2 + 2 



2 



where '^{z) is the imaginary part of z. We see that Qa > as x — 0. Thus we may rewrite (543) 
as 







k{z) = I -A[ — j ^i-^-^]nr nn{r,z)dr 



■X^l^^J i rn.(r,.)dr 



where in the last step we have replaced the integrand by the total mass concentration, Mq kg m~^. If 
the cloud particles are liquid droplets then the density p and index of refraction n are well known 
and Mq is referred to as the liquid water content. Assuming the extinction by cloud droplets 
dominates the volume extinction coefficient then we may integrate over the thickness of the cloud 
to determine the total optical depth. Referring to (545c) we have 



= (236) 
pA +2/ 

where Mq is the path-integrated water content in the cloud, called the liquid water path. It is 
notable that fe explicitly depends on the total condensed water mass, but not on the cloud droplet 
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size. In contrast, (95) shows that in visible wavelengths depends on both Mq and r. If the 
microwave optical depth of the cloud is known, through any independent measurement, then we 
may use it to obtain Mq by inverting (236). Fortunately fc{\) may be independently estimated 
through estimating the microwave brightness temperature of clouds. 

3.2 Anomalous Diffraction Theory 

The full Mie theory solutions are amenable to approximations in certain limits, the Rayleigh ap- 
proximation (§3.1) for example. Scattering and absorption by particles with indices of refraction 
near unity (n ^ 1) may be treated with anomalous diffraction theory (ADT) in the regime of large 
size parameters x ^ 1. ADT, described in van de Hulst [1957], estimates particle-radiation inter- 
action in this regime by considering the effect refraction has on the phase delay between scattered 
and diffracted waves. Consider incident plane waves of the form 

Ec(r,t) = Eoexp[i(k ■ r - wt)] 

= Eo exp[i(k ■ r — iujt + 0)] 

E(r,t) = g?(E,) 

= Eq cos(k ■ Y — ut + (p) 

= Eq^x cos(k ■ r — ut + (j)x)i + E^^y cos(k ■ r — ujt + 0y)j + i^o.z cos(k • r — ut +{^'Sjh 

where the phase angle cj) is an arbitrary angular offset. Of course plane waves which explicitly 
include (j) satisfy Maxwell's equations (474). is often omitted from plane wave solutions since, 
being a constant offset, it has no effect on the physical properties of the solution. We include in 
(237) because it will prove useful to examine the phase delay which accrues in plane waves due to 
particle scattering. 

After particle scattering, the electric field is comprised of two waves: the incident plane wave 
Eo which travels unrefracted through the particle (because n ~ 1), and the scattered wave. Eg 
which has been diffracted by the edge of the particle. 

E = Eq + E, (238) 

We assume that rays suffer no reflection or deviations at the particle boundaries since n ^ 1. Thus 
the particle changes Eg only in phase, not in amplitude. The change in phase may be determined 
by considering the geometric path length through the particle. CoUimated rays passing through the 
particle of radius r a distance a (the impaction parameter) from the center define chords of length 
2a sin ip from the entry point to the exit point. The angle ip lays between the radius to the point 
where the light ray enters the sphere, and the line segment (of length a) joining the center of the 
sphere to the chord the ray travels through the sphere. 

ip = cos~^ a/r (239) 

Thus the angle between the entry and exit points and the center of the particle is 2^. The phase 
lag A0 is determined by the phase of the electric field measured on a plane behind the particle but 
normal to the propagation of the incident wave. The electric field measured within the geometric 
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shadow of the particle has suffered a phase lag relative to the field outside the geometric shadow, 
which has not interacted with the particle. The path length traveled within the particle is a 

a = 2rsmilj (240) 

The phase lag A0 suffered during particle transit is a times the difference in propagation speeds 
times the spatial wavenumber k = 271/ A (12). We assume the particle is suspended in air so the 
propagation speeds differ by a factor of n — 1. 

A0 = a X (n — 1) X k 

= 2r sin ip x {n — 1) x k 

= psinip (241) 

where p is the phase delay suffered by a ray crossing the full diameter of the particle 

p = 2r(n- l)(27r/A) 

= 2x(n-l) (242) 

The phase delay (241) causes a radially varying reduction the electric field amplitude behind the 
particle, so that 

E, = E^e-'^^-l (243) 

We shall determine the extinction efficiency due to the particle by examining the scattering 
amplitude function S{Q) at the scattering angle = 0. The amplitude function 5(0) is the integral 
over the geometric shadow region of the ratio of the scattered wave to the incident wave. 

Following van de Hulst [1957], we employ polar coordinates within the geometric shadow region 
and begin by integrating over the impact parameter (239) 

a = r cos ip 

from the center of the particle to the edge. The element of area is a ring of radius a about the center 
of the particle. The elemental area is 2na da and area 

^^^^ = ^£(l-e-^^''^^)(27ra)da 

;i -e-'^^''^'^)ada (244) 



The first term in the integrand represents the incident wave. The exponential term which represents 
the scattered wave is more difficult to integrate. We now change integration variables from impact 
parameter to entry angle using (244), 

a = r cos ip 
da = —r sin ijj d?/' 
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The change of variables maps a G [0, r]tO'ip e [0, 7r/2]. Thus 

f3t/2 



3(0) 



1 2 2 

k r 



7t/2 



-ipsin^^lj^^ cos7A)(r sinV^)dV^ 



(245) 



where w = ip. The first term in the integrand in K{w) (245) is a simple trigonometric function. 
The second term may be integrated by parts with the formula j udv = uv — j vdu where 



u = 
dv = 

V = 
du = 

Using these standard techniques we obtain 

rn/2 



sin ijj 

„—w sin ?/) 



W 

COS Ip dip 



(246) 
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■ ■ 2- (247) 

With K{w) known we can finally determine the scattering function and the extinction efficiency in 
the anomalous diffraction limit 



^(0) 

Qe 



^$R[S(0)] 
A^[K{w)] 



For non-absorbing spheres, n is real so (248) becomes 

Qe 



A 4 , 

2 sinp H — ^-(l — cos p) 

p p^ 



(248) 
(249) 



Following a similar procedure \Liou, 2002] leads to the ADT approximation for absorption effi- 
ciency 



1) 



(250) 



where h = Axn-i. 
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3.3 Geometric Optics Approximation 

Nussenzveig [2003] estimate the fraction of absorption due to Mie resonance effects. First, they 
summarize previous studies [Nussenzveig, 1992; van de Hulst, 1957] which show Mie scattering 
in transparent spheres results in quasi-periodic resonances with approximate period 5x- For water- 
based liquid aerosols typical of the atmosphere, the period is of order unity in size parameter space. 
A more precise estimate is 

/X = - 1 (251) 

arctanu 

5x ~ (252) 

Nussenzveig [2003] estimate the absorption due to resonance effects based on complex angular 
momentum (CAM) theory. Let R be the fractional contribution of resonances to the total mean 
density of states. In the asymptotic limit as x ^ 1, 

R = ifi/n,f, forx<l (253) 

For H2O, Uj. = 1.33 so i? ~ 29%. Let ((5a,r) be the mean absorption efficiency due to resonances 
across the approximate resonance period Sx- Nussenzveig argues that, since each resonance adds 
to the absorption efficiency, i? is a good estimate of {Qas)/Qa., the fractional contribution of reso- 
nances to the overall absorption efficiency. 

A second estimate of {Qa,r)/Qa.,T can be obtained from the geometric optics approximation 
(GOA). A fundamental result of GOA [Bohren and Huffman, 1983] is that 

g 

ga,g = -n^A^-RhiX, fornix < 1 (254) 
Nussenzveig [2003] combines (254) with Wentzel-Kramers-Brillouin (WKB) theory to show 

(Qa,r) 3 arctan /i 



<5a,g 4 n3 - /i3 



2 

/i 



arctan /i 



1 



forx>l,ni<l (255) 



For TT-r = 1.33, (255) yields (Qa,r)/Qa,g ~ 16%, similar to the estimate based on R (253). 

Dave recommended a size parameter resolution for Mie integrations of Ax ~ 0.1. Markel 
[2002] shows this leads to an underestimate of the absorptance of soot in water spheres by about a 
factor of two. 



3.4 Single Scattered Intensity 

In clear sky conditions it is reasonable to assume that most of the majority photons measured 
by Downward looking satellite instruments measure upwelling radiation reflected by the surface 
plus any photons scattered by the atmosphere into the satellite viewing geometry. In clear sky 
conditions it is reasonable to approximate the solar radiance measured by a satellite as consisting 
entirely of singly-scattered photons. 

The upwelling intensity measured by a downward-looking instrument mounted, e.g., on an 
aircraft or satellite, is formally described by the exact solution for upwelling radiance presented 
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in (105). The full, multiple scattering source function with thermal emission S:^ (153) is difficult 
to obtain and we shall instead make the single-scattering approximation. In the single- scattering 
approximation, we replace S^{t, tl) by the single-scattering source function S'*(r, tl) (159). 

In the single- scattering approximation, the upwelling intensity leaving the surface consists 
solely of reflected photons (194) from the direct downwelling beam (150) 

i+{t\(i) = f ^i-{t*, n')p{u,-Ci',n) cos 9' dn' 

J — 

J — 

= /.oi^0e-^*/'^V(i.,-i7Q,f2) (256) 



3.5 Satellite Orbits 

Assume a satellite is in a circular orbit about Earth so that its Cartesian position r at radial distance 
r as a function of time t is 

r = r(sin wti + coswtj) (257) 

where u is the angular velocity. As its name implies, u measures the rapidity with which the 
angular coordinate changes (in radians per second, denoted s~^). Let r denote the orbital period 
of the satellite. We assume the orbital speed f of a satellite in a circular orbit is a constant. Circular 
geometry dictates that 



By definition 



T = 2%r/v (258) 



2n ■ 



= v/r (259) 

The acceleration a of the satellite is the second derivative of the position. 

V = ru {cos uti — sin utj) (260) 

a = — ra;^(sin wti + coswtj) (261) 

= -ioh (262) 

It is easy to show that v ■ r = because they are orthogonal. Thus the velocity of an object in a 
circular trajectory is tangential to the radial vector. In plane polar coordinates, r = rf so 

a = — rw^f (263) 



Moreover, (263) shows that the acceleration of an object in a circular trajectory is opposite in 
direction to the radial vector. 
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Since the orbit is circular, the magnitudes of r, v, and a are all constant with time and only 
their direction changes. This can be verified by computing the vector magnitudes, e.g.. 



a = |a| = v^a ■ a 

= rcj^ V sin^ cot + cos^ cot 

= ru^ (264) 

a is known as the centripetal acceleration. Substituting (259) into (264) we obtain 

a = r{v/rY 
= v'^/r 

a = f (265) 

r 

According to Newton's Second Law (??), the gravitational force F (??) holding the satellite 
in orbit must exactly balance (265). Let m and denote the mass of the satellite and of Earth, 
respectively. Then 

F = ma 

GmMff, . v'^ . 



— m — r 



— — — = m— 



r 



r 



^ (266) 

r 

We note that r is the distance from Earth's center of mass, not its surface, to the satellite. Hence 
satellite's speed increases slowly as its distance from Earth decreases. The orbital period r is 
determined by substituting (266) into (258) 



r = 2nr 



GMff 



Hence the satellite period squared is proportional to the cube of the orbital size. This relationship 
is known as Kepler's Law in honor of its discoverer, the astronomer Johannes Kepler. 



3.6 Aerosol Characterization 

The wavelength and size dependence of aerosol optical properties are generally well-predicted if 
the aerosol composition and shape are known. However, in situ aerosol composition and size are 
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often difficult or impossible to obtain in the field. Thus indirect estimates of particle size and 
composition are often derived from the measurements that are available. Many field sites and ex- 
periments are equipped to measure the spectrally resolved, column-integrated, aerosol extinction 
optical depth r^^^(\i) where i denotes the spectral channel. The following sections discuss tech- 
niques for inferring aerosol optical depth from surface measurements, and the subsequent uses of 
T^^'^{\i) to calibrate radiometry, and to infer information about aerosol size distribution and com- 
position. 

3.6.1 Measuring Aerosol Optical Depth 

Surface radiometry is incapable of directly measuring aerosol optical depth (AOD)^^ because 
aerosol is vertically distributed in the atmosphere. We now describe how AOD is inferred from 
surface measurements. This will enable us to qualify what is usually meant by the phrase "mea- 
sured optical depth" as it appears in the literature. Upward pointing radiometers may measure 
one or more parts of the radiance field. Downwelling total solar irradiance is traditionally 
measured by global solar pyranometers. The downwelling direct beam irradiance is most ac- 
curately measured by a normal incidence pyrheliometer (NIP). Downwelling diffuse irradiance F^ 
is measured by shaded pyranometers. Instruments that measure all three irradiances also exist. A 
class of instruments called shadowband radiometers do this by periodically blocking the detector 
from the direct solar beam using an occluding device called a shadowband. 

3.6.2 Aerosol Indirect Effects on Climate 

Greater than the direct radiative effect of aerosols on climate is the uncertainty are the effects on 
climate of changes in cloud due to aerosols. These changes are collectively known as Aerosol 
Indirect Effects (AIE). The first AIE is the droplet size effect or albedo effect described by Twomey 
[1977]. The second AIE is the so-called cloud lifetime effect or LWC effect. Other AIEs include 
the cloud glaciation effect. 

Dust may suppress precipitation [Rosenfeld et ah, 2001; Rudich et ah, 2002] including tropical 
storms and hurricanes [] or enhance precipitation [Levin and Ganor, 1996; Wurzler et ah, 2000]. 

3.6.3 Aerosol Effects on Snow and Ice Albedo 

Snow and ice always contain a number of particulate and bubble inclusions. The effect on re- 
flectance is non-negligible when the particles are substantially darker than the snow or ice. Warren 
and Wiscombe [1980] and Wiscombe and Warren [1980] show that Mie theory is an appropriate, 
but not exact, tool with which to model these effects. Light et al. [1998] present updated data on 
these effects. 

3.6.4 Angstrom Exponent 

The sensitivity of the particle extinction efficiency to wavelength i.e., \dQc{r, \)/d\\, generally 
increases with decreasing particle size. This sensitivity holds true any paritcle composition, and 

"AOD implicitly refers to aerosol extinction optical depth rj"^"^, the sum of aerosol absorption optical depth rf^"^ 
and aerosol scattering optical depth t^^^. 
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may be easily demonstrated using Mie codes. An empirical measure of this sensitivity is obtained 
by defining the Angstrom exponent ox Angstrom parameter a. 

ln[re(Ai)/re(A2)] 
ln(Ai/A2) 

Typical values of a are a > 2 for small carbonaceous aerosols (smoke) and urban pollution (sul- 
fate) to a ^ for large dust particles. 



4 Gaseous Absorption 
4.1 Line Shape 

Three fundamental processes interact to determine the observed shape of spectral lines. These 
are the natural line width, pressure broadening, and Doppler broadening. Natural line width, 
which arises from the Heisenberg uncertainty principle, is unimportant in planetary atmospheres. 
Pressure- and Doppler-broadening effects dominated line shapes here. However, the processes de- 
termining natural line widths and pressure broadening of lines are statistically identical, so these 
two processes have the same analytical form. In planetary atmospheres, Doppler broadening is 
intermediate in importance between natural line shape and pressure broadening or collision broad- 
ening. At pressures weaker than about 1 mb collisions become infrequent enough that the line 
shapes transition from Lorentzian to Doppler. In this transitional regime the line shape is a convo- 
lution of Lorentzian and Doppler shapes. This convolution is known as the Voigt line shape. 

4.1.1 Line Shape Factor 

The fundamental property of a discrete line transition is the spectral absorption cross section, 
a{h'), which measures the strength of the transition as a function of frequency (or wavelength, 
etc.). As discussed in Section 2.1.13, the units of a{u) depend on the application, and they must be 
consistent with the units of the absorber path, U . The absorber mass path is integral of the absorber 
amount u along the path. 



f/(Pl,P2) = 




(269) 



Thus, if a(z/) is expressed in units of m^ kg~^, then U and p should be expressed in units of kg m~^ 
and kgm~^, respectively. We adopt the convention that, when referrring to gaseous absorption, 
a{ij) is the number absorption coefficient and is expressed in m^molecule^^ Closely related is 
the volume absorption coefficient k.^ [m~^] 

Up) = a(z/)iVo (270) 

where Nq is the number concentration [molecule m^'^] of the species. 

The value of «(z/) depends on the absorber type. For a single spectral line, «(//) is obtained by 
converting the parameters in the hitran compilation [Rothman et al, 1998] to the local temper- 
ature and pressure. For a continuum absorption process, e.g., O3 Chappuis bands, a^v) is usually 
obtained from a standard table (which may include a parameterized temperature dependence). The 
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JPL compilation [DeMore et al, 1994] is the standard reference for absorption cross sections of 
photochemical species. 

Given these definitions, the optical depth due to absorption between points Pi and P2 is 

r^{.Pi,P2) = / a{iy)u{s)ds 



J Pi 

= a(z/) / u{s) ds 
Jpi 

= aiu)U{Pi,P2) (271) 

Comparing this to (82) we see that a{u) is the extinction due to absorption. 

The integrated strength of a transition is called the line strength or line intensity and is denoted 
hyS 

POO 

S= a(z/)dz/ (272) 



^0 

The dimensions of S are thus the dimensions of the absorption coefficient times the dimensions of 
frequency (or wavelength), e.g., m^ kg^^ Hz or m^ kg'^ m. The relation between a{u) and the line 
strength is called the line shape factor $(z/ — uq) 

a{iy) = S^{u - uo) (273a) 
S = a(z/)/$(i/ - z/o) (273b) 
$(zy - uo) = a{u)/S (273c) 

where z/q is the central frequency of the unperturbed transition. Thus $(z/ — z/q) relates the 
frequency-integrated absorption amount (i.e., line strength) to the specific absorption a distance 
— z/q from the line center. The shape factor for each transition (i.e., line) is normalized to unity 

<^{u-uo)du = 1 (274) 

The dimensions of $ are inverse frequency or wavelength, e.g., (s^^)^^, Hz^^, or m^^. 

In this section we do not need to specify the which set of self-consistent dimensions we choose 
for u, U, and a(z/). Nevertheless, for concreteness we list the three most common choices. First, 
we may employ dimensions relative to absorber number concentration 

u = n molecule m^^ 

U = N molecule m-2 (275) 

a{iy) = ks,{iy) m^ molecule"^ 

where ki,(u) is a number absorption coefficient. Alternatively, we may employ dimensions relative 
to absorber mass concentration 

u = m kgm~^ 

U = kgm-2 where 

a(z/) = ^p{l') m^kg"-*^ 
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iplu) is a mass absorption coefficient. (276) 

Finally, we may employ dimensions relative to absorber number 

u = m kgm~^ 

U = kgm-2 where 

a{iy) = ip{iy) m^kg"-*^ 

tpiu) is a mass absorption coefficient. (277) 

We prefer to work in units absorption per unit mass (277) because, at least in models, the mass of 
absorbers is more readily available than the number of absorbers. 

4.1.2 Natural Line Shape 

Radiative transitions from an upper to a lower state may occur spontaneously, in a process known 
as spontaneous emission or spontaneous decay. The time intervals t between spontaneous emission 
are described by a Poisson distribution 

pp{t) = -e-*/^^ (278) 

where is the mean lifetime of the excited state. 

The time-dependent Schroedinger equation predicts that the relation between the time intervals 
between such decays and the energy of the decays gives rise to a continuous, rather than discrete, 
profile of absorption. The resulting profile is called the natural line shape and is described by 

^nii^-^o) = -r. ^ (279) 

where is the natural line shape half width at half-maximum (HWHM). This is seen by noting that 
$n(z/ = z/q) = is the full maximum value of $n(^), while ^a{v = z/q + ^n) = (27Tan)~^- 

Natural broadening is one of two important broadening processes that are described by the 
Lorentz line shape, $l(z^ — t'o)- The Lorentz line shape is defined identically to (279) as 

*l(^-^o) = -r, ^ (280) 

The Lorentz profile (280) may also be written in terms of the full width at half-maximum 7l = 2an 

as 

*l(^-^o) = -7^ ^^Jl, ,..,1 (281) 

The HWHM of a Lorentz profile is often written in terms of a parameter Fl such that 

^^i^-^^) = ^( ^2^^"^^ /A^^2^ (282) 

n[{iy - voY + (FL/47r)^J 
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The full meaning of Fl is described below in §4. 1 .3. 

The Lorentz line shape describes the relative absorption as a function of distance from line 
center. A separate parameter, the line strength, S (272), defines the abolute absorption. Thus 
$l(z^ — vq) must be normalized such that its integral over all frequencies is unity. 
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We now make the change of variable 
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This change of variables maps v E [0, oo) to y E [—vo/an, oo). Since vq ^ 1, and the integrand 
is peaked near the origin, we will replace the lower limit of integration by — oo. 
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(283) 



Thus the Lorentz line shape (280) is correctly normalized. 

Some insight into the origin of natural broadening may be obtained from simple examination 
of the Heisenberg uncertainty principle which may be expressed as the fundamental uncertainty 
relating two conjugate coordinates such as position and momentum, or, in our case, time and 
energy: 

AtAE ~ h (284) 

When a molecule is in an excited state it has a probability of spontaneously decaying to the lower 
state. Its actual state is uncertain on the timescale At, during which it may be in a superposition of 
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allowed states. Corresponding to this uncertainty, then, is an uncertainty in energy AE = hAu so 
that 

AthAz/ ~ h/27l 

Au ~ (285) 
2TiAt 

If we identify At with Tg, the mean hfetime of the excited state, and Av with an, the natural half 
width, then 

«n = (286) 

2'KTe 

For spontaneous emission, Tg is the reciprocal of the Einstein A coefficient [Shu, 1991]. A typical 
value of Te is 10^^ s. This is an extremely long time relative to the typical mean time between 
(optical) collisions in a planetary atmosphere, order 10^^° s. In practice, therefore, the line shapes 
we are concerned with are dominated by coUisional effects. 

4.1.3 Pressure Broadening 

Collisions between molecules are the most important cause of line broadening in the lower atmos- 
phere. Statistically, collisions in a gas of molecules of uniform speeds take place at random time 
intervals about a mean value. As with spontaneous decay, the probability distribution of time in- 
tervals between collisions is therefore a Poisson distribution (278). In reality, molecular velocities 
are not uniform but obey the Maxwell distribution 

Pm{v) is the probability that a molecule of mass M at temperature T has a velocity in [f , f + df]. In 
all other respects, collision-induced line broadening, or pressure broadening of lines is equivalent 
to (279) and is thus governed by a Lorentz line shape (280). 

%{y-yo) = -r, (288) 

where the half width at half-maximum intensity due to pressure (collision) broadening, is called 
the pressure-broadened line width. Because Up charaterizes coUisional interactions, its value de- 
pends on the concentrations and masses of all the molecular species which are present. Thus its 
determination is somewhat involved. 

Molecular collisions occur in two senses, the kinetic and the optical. A collision might change 
the kinetic energy of one or both of the molecules, but this is not neccessary. Collisions which do 
alter the kinetic energy of one or both of the particpants are kinetic collisions. Collisions may also 
cause a radiative excitation or de-excitation of the radiative energy levels of either or both of the 
molecules involved. CoUisional de-excitation, for example, might release one vibrational quanta 
of energy without changing the kinetic energy of either of the molecules. Collisions which alter 
the rotational or vibrational quanta of one or both of the particpants are called optical collisions. 
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Rotational energy levels are, in general, much smaller than translational energy levels and thus 
optical collisions may occur more frequently than kinetic collisions. 

Many methods of determining the pressure broadened line shape have been proposed. A par- 
ticularly satisfying model is that of molecules radiating with constant intensity between elastic col- 
hsions that change the phase of the emitted wavetrain randomly. The radiation is a delta-function 
in frequency with harmonic time dependence 



fit) 



(289) 



uq is the center of the unbroadened line. We may restate (289) in frequency space by taking its 
Fourier transform 



Fiu) 



(290) 
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(291) 



The phase of the oscillation, and thus of f{t) is reset randomly by each collision. The waiting 
period between collisions determines the length of each wavetrain. The probability of a collision 
occuring between time t and t + dt follows a Poisson distribution 



3*/^^ dt 



(292) 



Let Tc and Tq denote the mean times between kinetic collisions and optical collisions, respec- 
tively. By definition the mean collision frequency Uco\ is the inverse of t^. The relationship between 
distance, rate, and time may be formulated to obtain 



Tr 



col 



(293) 



where is the thermal speed and A is the mean free path, both of which are known properties of 
the thermodynamic state of the atmosphere. First, let us recall that Maxwell-Boltzmann statistics 
prescribe a square-root dependence of v^a on T 



8RT 



8R*T 
nM 



(294) 



where Ai is the mean molecular weight of the gas, R* = A^i? is the universal gas constant, R is 
the specific gas constant, and T is the ambient temperature. Kinetic theory tells us [Seinfeld and 
Pandis, 1997, p. 457] that the mean free path Aaa of molecular species A in a gas of molecular 
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species A is inversely related to the total concentration of molecules and to their cross-sectional 
area for collisions, 7ro"^/4 

= 2, o (295) 
M 

MET 

R*T 
ny/2pj\fal 



(296) 



where'^ we have re-expressed Aaa in terms of p and T by using = pM / Ai and then applying the 
ideal gas law p = p/(RT) (??) to (295). There are many subtle assumptions embedded in (295) 
that should be clarified before proceeding: First is that Aaa is the mean free path of molecular 
species A in a gas of molecular species A, i.e., of a homogeneous gas of A. Often A is a trace 
species (e.g., CO2) in air, represented by species B. The mean free path of A in B incorporates 
properties of both A and B [Seinfeld and Pandis, 1991, p. 457]: 

The pressure of an inhomogeneous gaseous mixture of A and B is due, mainly, to the presence 
of the bulk medium which, for our purposes is air represented by species B. Although the tem- 
peratures T of both gases in a mixture are equal in thermodynamic equilibrium, their densities and 
partial pressures are unequal. Hence the expressions for pressure -broadened line shapes of A will 
contain molecular properties of A as well as bulk thermodynamics state (temperature, pressure) 
due mainly to B. 

For a homogeneous gas of A, substituting (294) and (296) into (293) yields 

8R*T li^pMol 
tlMI ^ R*T 
8R*T li^pMol 



4Ar, 



71 p 

R*MaVT 

AUai nr p 



Mk \ RaVT 



(297) 

We see the interesting result that ap decreases as T increases. This is because although thermal 
speed increases as y/T (294), the mean free path A between collisions increases linearly as T (296). 
The net result of increasing T is therefore to decrease collision frequency Ucoi, and ap (297). 



'^Unfortunately, terminology has placed in close proximity two symbols that are easy to confuse; A'a [#m is 
the number concentration of molecular species A, while J\f [#mol^^] is Avagadro's constant. 
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Table 5: Temperature Dependence of a/ 



Molecule 


Min 


Mean 


Median Max 


H2O 


0.28 


0.66 


0.97 


CO2 


0.49 


0.75 


0.78 


O3 


0.76 


0.76 


0.76 0.76 


N2O 


0.64 


0.78 


0.82 


CH4 


0.75 


0.75 


0.75 


O2 


0.63 


0.71 


0.74 


OH 


0.50 


0.66 


0.66 


SO2 


0.50 


0.60 


0.75 



"Sources: Rothman et al. [1998], Liou [1992] 



If we define 



ttp,0 = 



Po 



R/1 



(298) 



then 



p_ 
Po 



To 
T 



(299) 



where n = ^. The exponent n determining the temperature dependence in (299) is called the line 
width exponent. Two of the key line parameters measured in laboratory experiments are ttp o and 
n. Both are tabulated in databases such as hitran (§5.5.1). Table 5 shows the mean, median, and 
range of the temperature-dependent exponent n for the pressure-broadened half width Up of various 
optically active gases. For most optically active gases, 0.6 < n < 0.8. This range differs consid- 
erably from n = \ derived in (299), which is known as the classical value of n. This discrepancy 
arises becauses the optical cross section a of the molecules have a temperature dependence. 

Table 6 shows the mean, median, and range of the pressure-broadened half width Up of various 
optically active gases. 

Pure kinetic theory must be combined with the optical cross section a to determine ap^^ 



Olr. 



2kT\ / 1 1 



.1/2 



(300) 



^^(300) is taken from Goody and Yung [1989], p. 99. The equation appears to contain typos and should be viewed 
qualitatively until it is checked against a trustworthy reference on kinetic theory. For example, the quadratic depen- 
dence on the concentration of the perturber appears to be typo. Also the optical collision diameter a is missing. 
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Table 6: Pressure-Broadened Half Widths 



Molecule 


Min 


Mean Median 


Max 


H2O 


0.0077 


0.071 


0.11 


CO2 


0.055 


0.071 


0.095 


O3 


0.049 


0.069 


0.084 


N2O 


0.069 


0.075 


0.097 


CH4 


0.018 


0.054 


0.16 


O2 


0.028 


0.043 


0.060 


OH 


0.040 


0.044 


0.095 


SO2 


0.10 


0.11 


0.15 



"Units are cm ^ atm ^ at 296 K. Sources: Rothman et al. [1998] 



where the summation is over all perturber species i (i.e., O2, N2, . . . ) and M is the mass of the 
absorber. Thus ap depends quadratically (fxm: typo, must be linear) upon the number density 
of collision partners and linearly upon the velocity of the molecules (note the factor of \/VT). A 
scaling approximation is therefore frequently used to extrapolate ap{T) based on some measured 
or tabulated value a;p(To) 

(rp \ 1/2 
^ I (301) 
To/ 

Atmospheric pressure changes by three orders of magnitude, from 1-1000 mb, from 50 km to 
the surface. Temperature changes by only a factor of 2 over the same altitude range. Since line 
widths (301) depend linearly on p, but only on the square root of T, pressure variation dominates 
vertical changes in ap. 

The mean time between optical collisions Tq is then related to Up by 

1 



2%ap 



ap 



o 



«p = ^ (302) 

This is exactly analogous to the relation between Tg and (286) because both processes are 
described by the Lorentz line shape. 

When natural and collision broadening are considered simultaneously, it can be shown that the 
combined line shape is also a Lorentzian with HWHM rL/47r (282). The parameter Fl is related 
to both the natural line width and the mean frequency of collisions Ucoi = (302) by 



(303) 



4 GASEOUS ABSORPTION 



79 



These relations may be proved by examining the power spectrum of a sinusoidal electric field 
which is randomly interrupted Vcoi times per second. In practice in the lower atmosphere ap >> 
so the Lorentz line width is nearly equivalent to the pressure -broadened line width. Therefore, 
Lorentz broadening and the Lorentz line width will hereafter refer to the Lorentzian line shape 
(280) due to the convolution of natural and collision broadening. With reference to (303) and 
(302) 

dh = ftn + i^coi/27r 

QiL = On + Ctp (304) 

When discussing a given transition, it is convenient to translate the origin of the frequency axis 
to the line center. Rather than defining a new variable to do this, we continue to use u. Thus we 
replace u — u^hy v alone. The intended meaning of v should be clear from the context. Using this 
notation (280) becomes 

+ al) 

The absorption cross-section (273) for Lorentzian lines may now be written 

n{u^ + al) 

4.1.4 Doppler Broadening 

Brownian motion causes molecules to constantly change their velocity relative to the photons 
which may interact with them. While the frequencies of resonant absorption are constant in the 
frame of motion of the molecule, this random motion broadens the range of resonant frequencies 
in the frame of a stationary observer. This form of line broadening is called Doppler broadening. 
The probability that the molecule and a stationary reference frame (an "observer") have relative 
velocity in [v,v + dv] is given by the Maxwell distribution (287) 



f M ( Mv' 

pyi{v)dv = exp --zTT^ dv (307) 



1/2 

\2KkT J \ 2kT 

We convert this Maxwellian PDF directly to the Doppler line shape PDF $d (z^ — t'o) by noting that 
each value of v describes a shift in line center /S.v = u — uq. For f /c ^ 1, the relation between v 
and the Doppler shift Az/ is 

Az/ = z/Qf/c 
V = cAu/uq 

dz/ = dAu = — dv 
c 

dv = — du 

^0 
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When we re-express Pm(^) (307) in terms of Au by imposing the condition 



(308) 



This relation between pm{v) and puiAu) is like that between Bx(\) and By{u) (42). Inserting 
(308) into (307) leads to 



Pm[v{v)] dv 



1/2 



\2t±T ) 
V27rkT ) 



exp 



Mc^Au'^\ dv 



dz/ 



2z/2kT 



du 



c 



1/2 



exp 



2z/2kT 
Mc2 



di/ 



(309) 



The RHS of (309) is now strictly a function of u. We name this function the Doppler line shape, 
^'^{y) and define the Doppler width ckd as the square root of the portion of the exponential in 
parentheses 



c 



so that (309) becomes 



$d(z^) du 



(- 



M \ 



1/2 



/ M 



1 c 
1 



2kT 



2kT 

~M 

exp 

1/2 



(310) 



Au" 

X.2 



a 



dv 



exp 



Av' 



du 



exp 



dz/ 



(311) 



The normalization of (311) follows from the fact that it is simply a re-expression of the Maxwell 
distribution (307), which is already known to be normalized. The full maximum value of $d is 
$d(z/ = z/q) = (ao-v^)"^, and ^-q^u = z/q + ckd) = {eaY^y/iiY^ , so that c^d is the half width 
at e"^ of the maximum of $d- Care should be taken not to confuse ckd with a half-width at half- 
maximum. The half width at half-maximum of is 



(312) 



It is more appropriate to compare q;l (298) to ao (312), since they are identical measures of line 
width, than to ao (310). 

We recognize that (31 1) is a form of Gaussian distribution, although it is not expressed in the 
canonical form of a Gaussian 



Pg{v) dv 



1 



(TV 271 



exp 



aV2 



dz/ 



(313) 



in which a represents the standard deviation and z/q the mean value. Comparing (313) to (31 1) we 
see that the standard deviation of the Doppler line shape ctd = od/ \/2- Thus Doppler broadening 
alone allows absorptions to occur within ctD / and -\/2«d of line center with efficiencies (relative 
to line center) of 1 — e~^ = 0.683, and 1 — e"^ = 0.954, respectively. 
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4.1.5 Voigt Line Shape 

Line broadening processes occur simultaneously in nature. In particular, Collision broadening 
occurs in tandem with Doppler broadening. If we assume these two processes are independent but 
occur simulateously, then the net line shape of the total process will be the collision-broadened line 
shape, shifted by the Doppler shift (308), and averaged over the Maxwell distribution (307). The 
resulting line shape is call the Voigt profile, $v(^)- 

r f M V'^ ( Mv^\ aL 

= / — — exp -— — dv (314) 

\2nkT J V 2kr J n[{iy - uq + vqv/cY + a£j 

-TTTTTT -FT \ TT^^ ^^17 (315) 



2nkT V 2kT y n[{u - z/q + I'ov/cy + "i 

Analytic approximations to and asymptotic behavior of $v(i^ — i^o) are discussed in Goody and 
Yung [1989]; L?oM [1992]. 

With appropriate definitions [e.g., Liou, 1992, p. 30] 

t = yin2(z/ - z/')/ttD 
\/ln2(z/ — h'o)/aB 



X 



n2aD 



y 



y/ln 2ap/aD (316) 



we may re-express (315) as 



$v(i^-^^o) = —\—K{x,y) (317) 

V TT 



where the complex error function K{x, y) is defined 

Efficient algorithms for evaluating (a;, y) have been developed [//M/e? a/., 191 S; Humbcek, 1982; 
Kuntz, 1997] to reduce computational expense in detailed line-by-line calculations. 



5 Molecular Absorption 

We now attempt to develop an introductory understanding of the location of absorption lines in 
gases. The science of atomic and molecular spectroscopy is extremely detailed and requires a 
quantum mechanical treatment for a satisfactory explanation of all phenomenon. Nevertheless, 
important insights as to the spectral location, temperature dependence, and relative population of 
radiative energy levels may be gained by resorting to a semi-classical treatment of the radiation 
field. In essence, our first task is to characterize the distribution of observed molecular transitions 
in the atmosphere, i.e., to characterize AE from Planck's relation 



AE = hi/ 



(319) 
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To begin, we enumerate the contributions to a molecule's total energy E^. The total molecular 

energy comprises the tranlational kinetic energy ngrtrn, the electronic energy E^, the vibrational 
energy Ey, the rotational energy Ey-, and the nuclear energy E^. 

Et = Et + E^ + Ey + E, + Er, (320) 

These energy components have been listed in order of decreasing magnitude. We are mainly 
interested in describing vibrational and rotational energy transitions in the next sections. 



5.1 Mechanical Analogues 

In the quantum mechanical view, radiative absorption occurs when an incident photon of energy 
level hu encounters a molecule with energy Ei. If Ei+hu is "close enough" to an available energy 
state E2 of the molecule, then absorption may occur. If absorption does occur, the energy of the 
photon is transferred to to the corresponding energetic mode or modes of the molecule. These 
molecular modes, e.g., spin, vibration, rotation, take their names from mechanical systems which 
form their classical analogues. To understand the distribution of energy levels in these modes it is 
thus useful to describe classical mechanical analogues which are known to behave like molecular 
systems. 

Molecules may be thought of as rotating structures. When the molecular structure is fixed, i.e., 
the separation between the atoms does not change, we call the molecule obeys the rigid rotator 
model. Diatomic molecules, for example, may be visualized as barbell-shaped, with an atom at 
each end, separated by a massless but rigid rod. The separation of the atoms means their charges 
are separated. If the center of mass between the atoms does not coincide with the center of charge, 
then the molecule has a permanent dipole moment. Molecules comprised of a single element, e.g., 
N2 or O2, are called homonuclear. Homonuclear molecules are perfectly symmetric rigid rotators, 
and thus have no permanent dipole moment. 

To continue with the quantum mechanical behavior of the rigid rotator molecule, we must 
determine the available energy states of the system. Let the atoms of masses Mi and M2 be 
separated by a distance r. Then the distances of each from the center of mass of the system are 

Mir M 

ri — — — = — — r (321) 



Ml + M2 M2 
M2r _ M 
Mi + M2~ Ml 



M2r M 

r2 = — — = —r (322) 



where the reduced mass M of the system is 

* - (323) 

Ml + M2 

The moment of inertia I of the system about the center of mass is the total mass- weighted mean 
square distance from the axis of rotation 



I = Miri+M2r2 



(324) 
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The angular momentum L of the system is the total mass-weighted product of the distance from 
the axis of rotation times the velocity. Denoting the angular frequency of the rotation by u, the 
linear velocity of the atoms are uri and ujr2- 

L = MiTiVi + M2r2V2 
= Miurl + Mswr^ 
= uj{Mirl + M2rl) 

The classical energy of a rigidly rotating dumbbell may be expressed in terms of / and L 

= ^ (325) 

The quantum mechanical angular momentum operator, C, quantizes the angular momentum 
component of the energy Hamiltonian. The eigenvalues of the square of the angular momentum 
operator, are 

£2^ = n^J{J+l)tlj (326) 

where ^ is the wavefunction and J, an integer, is the rotational quantum number of the system. 
Substituting (326) into (325) we obtain the quantum mechanical energy of a rigid rotator 

1 c 

J(J + 1) 



47l2 21 c 

-^hcJ(J + l) 

hcB^J{J+l) (327) 



where 



i?v is the rotational constant of the species. is purely a function of the atomic masses and 
geometry of the species through I (324). The v subscript''* indicates that B^ is a function of 
vibrational state. 

Spectroscopists have adopted a variety of equivalent but usually unintuitive notations to de- 
scribe the lower (lower energy) and upper (higher energy) states of transitions. The most common 
convention is that lower and upper state quantum numbers are superscripted by a double prime 
""" and a single prime respectively. Thus lower and upper vibrational and rotational quantum 
states are denoted by {v", J") and (f ', J'), respectively. Examples of equivalent representations of 



By is not to be confused with the Planck function, B^. 
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emission include 



{v',J') -> {v",r)+hiy 
{v",J")+hu ^ {v',J') (329) 

Due to destructive interference patterns between quantum mechanical wavefunctions, not all 
conceivable vibrotational transitions are allowed transitions. Rotational transitions are subject to 
the selection rule that 



A J = J'±J" = ±1 (330) 

Transitions which break selection rules like (330) are called forbidden transitions. 

Applying this rule to (327) we obtain the energy of photons released in purely rotational emis- 
sion 

AE, = hc5^(J + l)(J + 2)-hc5^J(J+l) 
hu = hc5v[j2 + 3J + 2- (J2 + J)] 
u = cSv(2J + 2) 

= 2c5v(J+l) (331) 

Thus lines in pure rotational bands are spaced linearly in frequency. Linear, symmetric molecules 
such as N2, O2, and CO2 have no permanent dipole moments in their ground vibrational states 
because the center of mass coincides with the center of charge. These molecules therefore have no 
rotational bands in their ground vibrational state. 

The factor of c on the RHS of (33 1 ) disappears if we work in wavenumber rather than frequency 
units. Using i> = u/c results in 

z> = AE, = 2B^{.J + 1) (332) 

This direct relation between wavenumber and quantum number is one reason spectroscopists prefer 
to work in wavenumber space. 



5.1.1 Vibrational Transitions 

A molecule composed of atoms has K vibrational modes where 

_ { 3N — 3 : Non-linear molecules f333) 
1^ 3A^ — 2 : Linear molecules 

Often modes are degenerate, and so have fewer distinct vibrational quantum numbers than indi- 
cated by (333). As shown in Table 7, linear molecules include all diatomic molecules (CO, NO, 
OH), as well as CO2 and N2O. Non-linear molecules include H2O, O3, CH4, and CFCI3. 
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Table 7: Mechanical Analogues for Radiatively Important Atmospheric Gases" 

Moments of Inertia Class Members 
/^ = 0,/b = /c;^0 "Linear CO2, N2O, CO, NO, HF, HCl, HBr, HI, OCS 

^ 0, Is = Jc 7^ Symmetric top CFCI3, NH3, CH3CI, CaHe, SFg 
Ia = Ib = Ic Spherical top CII4 
I^^I^^ Ic Asymmetric top H2O, O3, SO2, NO2, HNO3 

"Based on Goody and Yung [1989], p. 81, Rothman et al. [1998] 
^Diatomic molecules and linear polyatomic molecules with integer J 

The equation of motion for the position r of a classical oscillator of mass M oscillating with 
restoring force j3 is 

M— = -/3r 



r 



A 

M 




(334) 

where (3 is the constant of the linear restoring force, known as the spring constant or the bond 
strength. Since (3 is the restoring force divided by the displacement distance, its SI units are N m^^ 
(Newtons per meter) and its COS units are dyncm^^ (dynes per centimeter). The spring constant 
(3 depends only on the total particle mass, not on the mass distribution within the particle. The 
solution to a two particle sysytem where particle equilibrium separation is r is identical to (334) 
with the mass M replaced by the reduced mass M (323) 




M1M2 

The bond strength (3 depends on the exact arrangement of nuclei within a molecule. For diatomic 
molecules (3 approximately fits the following progression: Single bond, double bond, and triple 
bond diatomic molecules have (3 of approximately 500, 1000, and 1500Nm~^(5 x 10^, 10 x 10^, 
and 15 x lO^dyncm^^) respectively. 

Let us consider the aborption spectrum due to stretching of the carbon monoxide molecule. 
CO is approximately a triply bonded diatomic molecule with a bond strength of about 1900 N m^. 
We will use (335) to predict the fundamental frequency vqq of the absorption spectrum. We shall 
include the details of the calculation to demonstrate the equivalence of the results obtained using 
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SI and CGS units. 




1900 Nm-i X (12.0 x 10-3 + 16.0 x 10-3)kgmori 
12.0 X lO-^kgmor^ X 16.0 x lO-^kgmor^ 



53.2 Nm-i / ,kgm 1 A/" 

- 0.160a / 2.77 X 105^:^ X - X 



192 X 10-6kgmori V m kg 



6 02 X 10^3 

z/co ~ 84.2W^x^ ^ 65.3 X 10^2 s-i ^ g5_3THz 

c 3.00 X lO^ms-^ R 

Ago = ^— ——— ^ = 4.59x10-^ m = 4.59 ^im 

z/co 65.3 X 10^^ s ^ 



65.3 X 10^^s~ 



z>co = ~ ] 5 T ~ 2182 cm-^ (336) 

100c 100 cmm-i X 3.00 X 108 ms-i ^ 

The conversion of uco to z>co and Aco is straightforward, and shows that the stretching mode of 
the CO molecule produces an absorption spectrum in the infrared. For pedagogical reasons, let us 
repeat the derivation of z/co using CGS units to demonstrate that neither system is clearly superior 
in terms of its simplicity or ease of manipulation. 




X lO^dyncm-i x (12.0 + 16.0) gmor^ 
12.0 gmor^ X 16.0 gmor^ 



192gmol W s cm g 

and the rest of the computation is identical to (336). It appears that CGS units, historically the 
system of choice for radiative transfer, are slightly more simple to work with than SI units because 
of the CGS-centric definition of atomic weights. However, in this era of automated computations, 
the slight advantage of CGS does not appear to outweigh the usefulness of consistently using SI 
units (336). All other things being equal, we recommend using SI units wherever practical in 
radiative transfer, a convention we follow in this text. 

The vibrational energy of a quantum harmonic oscillator is 

K 

= ^hz/fc(i;fc + i) (337) 

k=l 

where z/^ is the frequency of the kth mode, hz/^ is the vibrational constant of the mode, and Vk 
is the vibrational quantum number of the mode. Like other quantum numbers, v is an integer 
V G [0, 1, 2, . . .]. The term of ^ in (337) is called the zero-point energy. The vibrational ground 
state of a molecule is reached when = for all K vibrational modes. The linear relation 
between Vk and (337) is an approximation. In reality anharmonicities cause oscillations to 
deviate from linear behavior. This causes vibrational energy levels to become more closely spaced 
as Vk increases. 

Allowed transitions must match applicable selection rules. The selection rule for vibrational 
transitions is 



Aw = ± i;" = ±1 



(338) 
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The fundamental transition refers to transitions between the first excited state (v = 1) and the 
vibrational ground state (v = 0). Using Av = ±1 in (337) leads to 

AE^ = hiyk{v + l + l)-huk{v + l) 
= h[i^kiv + I) - Ukiv + I)] 

= hz/fc (339) 

For simpler molecules such as diatomic molecules, the relation between the type of vibrational 
mode k (e.g., stretching, bending) and the fundamental energy of the mode (339) may depend on 
only a few parameters. 



5.1.2 Isotopic Lines 

The quantum mechanical analogue of (334) is 




(340) 

where now (3 is the force constant of vibrational mode k which depends on the distribution of 
charge in the molecule but is completely independent of the mass distribution, and thus independent 
of isotope. Two distinct isotopes i and j of the same species will thus have differing equilibrium 
frequencies Ui and Uj for the same vibrational mode: 



- = \h7 (341) 
For instance, the frequency shift for ^^C^^O relative to the 2140 cm^^ band of ^^C^^O is 47.7 cm^^. 



5.1.3 Combination Bands 

Transition selection rules for many important molecules (including all diatomic molecules), require 
that Av 7^ and A J 7^ 0, i.e., vibrational transitions must occur simultaneously with rotational 
transitions. Thus simultaneous transitions of both vibrational and rotational states are the norm, not 
the exception, in the atmosphere. These simultaneous transitions form what are called combination 
bands or vibration-rotation bands. The energy of such transitions is obtained by subtracting the 
lower energy state from the upper energy state 

PR = Pk + B'J'{J' + l)-B';j"{J" + (342) 

where Uk is the energy of the pure vibrational transition (337). We have not cancelled like terms, 
as in the idealized case of (331), since B'^ ^ B'^ in combination bands. Substituting A J = ±1 into 
(342) leads to the energy spacings between the R-branch transitions and the P-branch transitions, 
respectively. 



i>R = h + 2B'^ + i3B'^-B':)J+iB'^-B':)j' 
up = h-iB'^ + B':)J+iB'^-B':)j' 



(343) 
(344) 
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There may be a number of absorption bands in addition to the fundamental bands of a molecule. 
These include overtone bands, combination bands, vibration-rotation bands, and harmonic cou- 
pling bands. Overtone bands occur at frequencies which are multiples of the fundamental fre- 
quencies. Combination bands are due to the interaction of two fundamental vibration bands. The 
combined frequency . . . Vibration-rotation bands have already been discussed. Harmonic coupling 
bands occur when interactions among closely spaced oscillation frequencies produces distinct, un- 
expected bands. This is relatively uncommon. 



5.2 Partition Functions 



The statistical weight g measures the number of available states in a given quantum configuration. 
Quantum theory teaches us that rotational energy levels are 2 J + 1 degenerate owing to 2 J + 1 
indistinguishable orientations of the component of angular momentum in a fixed direction in space. 
The statistical weight for rotational levels is therefore 



9r 



2J+ 1 



(345) 



In LTE, rotational states are populated according to Boltzmann statistics, so that the probability 
of occupancy is proportional to e^^"/'^^. At temperature T, the ratio of molecules in state J' to 
those in state J" is therefore 



N{J') 
N{J"] 



, Q-Er(J')/kT 



9 J" e 



-E,{J")/\T 



2J' + l 
2J" + 1 



exp 



[J'{J' + l)-J"{J" + l)] 



(346) 



where now J' and J" refer to any values of J, not just those satisfying selection rules. To examine 
the fractional abundance of molecules in state J relative to those in all other rotational states, we 
must sum (346) over all J" while holding J' = J fixed. This leads to 



N 



2J+1 



Qr 



exp 



hcB 



kT 



-J{J + l) 



(347) 



where we have defined the rotational partition function Qr as the total, probability- weighted num- 
ber of available rotation states in the system 



Q, = ^(2J+l)exp 



j=i 



J(J+1) 



(348) 



For atmospheric cases of interest, kT ^ hcB^ so the exponential term becomes vanishingly small. 
To express and evaluate (348) as an integral, we make a change of variables 



X 



dx 

(2J+ l)dJ 



hc5, 
kT 



J(J + 1) 
(2J+ l)dJ 



-dx 



(349) 
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This maps J G [0, oo) to x G [0, oo). Thus x is considered a continuous function of J. Substituting 
(349) into (348) results in 



Qr 



(2 J + 1) exp 

' kT 

-e dx 



' kT 



dJ 





kT 
kT 
kT 



hcB^ 



[-o-(-i: 



Substituting (350) into (347) we obtain 

N{J) hc5v(2J + l) 



N 



kT 



exp 



hcB^ 



J{J 



(350) 



(351) 



Examination of (351) shows that N{J)/N is small for J —i- and for J — )• oo and maximal in 
between. The derivation of (351) neglected any vibrational-dependence of B^. 

The vibrational partition function is defined by a procedure analogous to (346)-(350). 
Contrary to rotational states (345), vibrational states all have equal statistical weights 



9v = 9o 

In LTE, Boltzmann's Law, (352) and (337) lead to 

N{v') 
N{v") 



(352) 



9v' 


Q-E,{v')/VT 








e-^v{f")/kT 








— exp<^ - 
9o I 










kT 




kT 


exp 


kT^ 


v' - 


v") 







2 



(353) 



Notice that the statistical weight and the zero-point energy contribution of hz/o/2 factor out of the 
relative abundance of molecules in a given vibrational level (353). Nevertheless, the zero-point 
energy does contribute to the total internal energy of the system and so is properly included in the 
definition of the vibrational partition function 



X^exp 



i;=0 



hz^ 
'kT 



oo 

v=0 
^-huo/2kT 

I — Q-huo/kT 



(354) 
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where in the last step we used the solution for an infinite geometric series whose ratio between 
terms r < 1, 

oo 

Vaor" = (355) 
1 — r 



The fractional abundance of molecules in vibrational state v is thus 



N(v) 1 



kT 



1 _ p-hi/o/kT 

^ ^ -hi/o/2kTp!;hi/o/kT 

g-hi/o/2kT 

= (1 _ e-h^o/kr)e-''^'^«/'^^ (356) 

Vibrational levels are so widely spaced that most molecules in the lower atmosphere are in = 
or t> = 1 states. 

It is common to approximate the total internal partition function of atmospheric transitions Q 
as the product of the vibrational and the rotational partition functions. Using (354) and (350) we 
obtain 

1 _ e-h,yo/kT hc5^ 

In the high temperature limit where kT ^ hz/g then we may use the behavior ^ 1 — x so this 
becomes 

1 - 1 + ^ hcB, 

^ 2kT 

^ hcB, 
2kT-hz/o kT 
2hi/o hcSv 

5.3 Dipole Radiation 

The probability Vif of absorption of a photon leading to a change in molecular state from state i to 
state / is 

2 

dfi (357) 



-P^f = ^(j^)E/K„(k)|(^,|e^'^-e.(k).p|^,)P 



,ik-r 



fi 



The matrix element (V^/|e ea(k) ■ p|V^i) determines the absorption probability. The dipole ap- 
proximation retains only the first term in the expansion 

e^'^-" = l + ik-r + --- 

so that the matrix element simplifies to 

(7/>;|e^'^-^e,(k)-p|^,) = e„-(^;|p|^,) (358) 

This matrix element is related to the line strength S as follows [Rothman et al, 1998, p. 709].... 
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n2 




ni + hz/Q 


ni + hi^o 




n2 


n2 + hi/Q 


B21^ 


111 + 2hz/o 


n2 + M + KE 




111 + M + KE' + lii/( 


111 + M + KE 


Ci^ 


112 + M + KE' 



5.4 Two Level Atom 

Einstein created an elegant paradigm for analyzing the interaction of matter and radiation. His 
tool is called the two level atom. It consists of an ensemble of molecules with two discrete energy 
levels, El and E2 > Ei. The energy levels have statistical weights gi and g2, respectively. There 
are rii molecules per unit volume in the state with Ei, and n2 in the state with E2. The frequency 
of transition between the two levels is \iuq = E2 — Ei. 

The five processes by which an idealized two-level molecule may change state are 

(359) 
(360) 
(361) 
(362) 
(363) 

where M is a collision partner and KE and KE' represent the kinetic energy of the two molecule 
system before and after the collision, respectively. The most frequent collision partners are, natu- 
rally, nitrogen and oxygen molecules. 

The Einstein coefficient y42i s"^ is the transition probability per unit time for spontaneous emis- 
sion. Spontaneous emission, as its name implies, requires no external stimulus. Thus decay of 
excited molecules occurs even in the absence of a radiation field. Each emission reduces the popu- 
lation of excited molecules and increases the population of ground state molecules by one. Hence 

^ = ^1^2 (364) 
at 

The Einstein coefficient B12 determines the rate of radiative absorption. B12 is also called 
the Einstein coefficient for stimulated absorption. More specifically, -B12 is the proportionality 
constant between the mean intensity of the radiation field and the probability per unit time of 
an absorption occuring. The transition probability per unit time for radiative absorption is -812/1/- 
B12 has units of [s^^(Wm^^ sr^^ Hz^^)^^]. The rate of radiative absorptions per unit time per unit 
volume is 

dn2 r f , ,h 



dt 



rii [ [ a(z/)— ^ di7 dz/ 
Jo Jn hz/ 



47lni / a(z/)-^dz/ (365) 
hu 



where a;(z/) refers to absorption cross-section of the single transition available in the two-level 
atom. Comparing the final two expressions we see that 

Anaju) 

^12 = 1 ^/ N (366) 
hz/$(z/) 



B12 is intimately related to the line strength S of the molecular transition [Rothman et ah, 1998, 
p. 709]. In particular, B12 may be expressed in terms of the weighted transition-moment squared. 
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The Einstein coefficient B21 determines the rate of stimulated emission or induced emission. 
B2ilu is transition probability per unit time for stimulated emission. B21 has the same dimensions 
as B12. 

^ = n2 r [ a%u)^dndu 
dt Jo Jn hz/ 

= Ann2 / a\v)-^dv (367) 
Jo 

where a^{v) is the absorption cross-section for stimulated emission. 

Gaseous absorption cross-sections a{u) are extremely narrow and so the integrals over fre- 
quency domain in (365) and (367) are determined by a small range of frequencies about z/q, the 
line center frequency. Exact descriptions of line broadening about z/q are discussed in §4.1.2- 
§4.1.5. It is instructive to ignore the details of the finite shape of fine transitions for now and to 
assume the line absorption cross section behaves as a delta-function with integrated cross section 

a{v) ai2(5(z/ - z/o) (368) 
ai2 = (369) 

where rrie is the electron rest mass, e is the electron charge, and / is the oscillator strength of the 
transition. With these definitions. 

In thermodynamic equilibrium the number of transitions from state 1 to state 2 per unit time per 
unit volume equals the number of transitions from state 2 to state 1 per unit time per unit volume. 
Combining (364), (365), and (367) we obtain 

drii dn2 



dt dt 

niBi2lu = ^12^21 + n2B2iIu (370) 



5.5 Line Strengths 

We repeat the definition of line strength (272) for convenience 



00 



S= / a{v)du (371) 
Jo 

In practice, (371) is either measured in the laboratory or predicted from ab initio methods. 



5.5.1 HITRAN 

The HITRAN ("high-resolution transmission") database provides the parameters required to com- 
puted absorption coefficients for all atmospheric transitions of interest. The HITRAN database is 
defined by Rothman et al. [1998], and its usage is defined in their Appendix A. 

The idiosyncratic units of line parameters can be confusing. Most HITRAN data are provided in 
CGS units. For reference. Table 8 presents the equivalence between the symbols used in Rothman 
et al. [1998], their Appendix A, and our symbols. 
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Table 8: HITRAN database' 



This Text 


HITRAN 


Description 






Line strength 




7air 


Pressure-broadened HWHM 




7self 


Self-broadened HWHM 






Transition frequency 


i>i' 




Lower state energy 


n 


n 


Exponent in temperature-dependence of pressure -broadened halfwidth 


6 


5 


Air-broadened pressure shift of transition frequency z>o 


V 


V 


Frequency in wavenumbers 



"Equivalance between symbols employed in this work and those used by Rothman et al. [1998], Appendix A. 
^HITRAN data are tabulated at a temperature Tq and pressure po of 296 K and 101325.0 Pa, respectively. 



The line strength Si^ (272) is tabulated in HITRAN in CGS units of wavenumbers times cen- 
timeter squared per molecule, cm^^ cm^ molecule^ ^, often written cmatm^^ Line strengths for 
weak lines can be less than 10^^® cm^^ cm^ molecule^^, i.e., unrepresentable as IEEE single pre- 
cision (4 byte) floating point numbers. Thus it is more robust to work with Sifl units of per mole 
(rather than per molecule) in applications where using IEEE double precision (8 byte) floating 
point numbers is not an option either due to memory limitations or computational overhead. Fol- 
lowing are the conversions from the tabulated values of Si^ in cm^ cm^ molecule^ to more use- 
ful units, including cm~^ m^ molecule^ \ cm^^ m^ mol^^, and to fully SI units of Hzm^ mol~^ and 
mm^ mol^^. 

S'i,o[cm~^ m^ molecule"^] = Sifi[cm~^ cm^ molecule-i] x 10"^ (372a) 
Sifl[cm-^ m2 mol-i] = Sifl[cm-^ cm^ molecule-^] x 10"^ x (372b) 
S'i,o[Hzm2 mol-i] = Sifl[cm-^ cm^ molecule-^] x 10"^ x 100c x Af (372c) 
5i,o[mm2 mol-i] = Si,o[cm'^ cm^ molecule-^] x 10^^ x lOOA^ x Af (372d) 

where J\f is Avagadro 's number. Equation (372d) is uncertain. 

The energy of the lower state of each transition, E'l, is required to determine the relative pop- 
ulation of that state available for transitions. HITRAN supplies this energy in wavenumber units in 
the tabulated parameter v'- which is simply related to E'l by 

E'l = hcD'l (373) 

Of course, consistent with HITRAN and spectroscopic conventions, z/f is archived in CGS, cm^^. 

All HITRAN data have been scaled from the conditions of the experiment to a temperature and 
pressure of 296 K and 1013.25 mb, respectively. Thus line strengths must typically be scaled from 
HiTRAN-standard conditions, S'j_o(To), to the atmospheric conditions of interest, Si{T). 

SiiT) = 5i,o(To) ^^^^ g-hci>,/kTo 1 _ e-hc^../kTo ^^^"^^ 
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where Q is the total internal partition function. The second, third and fourth factors on the RHS 
of (374) are the temperature-dependent scalings of, respectively, the total partition function, the 
Boltzmann factor for the lower state, and the stimulated emission factor. The scaling presented in 
(374) is general and contains no approximation. 

The classical approximation [Rogers and Yau, 1989; Thomas and Stamnes, 1999] assumes that 
Q is the product of the rotational and the vibrational partition functions 

Q{T) ^ Q,{T)Q,{T) (375) 

In this approximation, (354) and Qj. (350) are treated as independent components of Q. This al- 
lows Q to be computed as discussed in §5.2. Most high precision calculation improve this approach 
by using more exact methods described in Gamache et al. [1990] and Gamache et al. [2000]. 
For many applications, a simplified version of (374) is employed 



SiiT) = Sifl(To) ( ^ ) exp 
^i,o(To) ( ^ ) exp 



k 











To. 


)] 


1 













(376) 



where m is an empirical parameter of order unity and The parameterization (376) incorporates 
the temperature dependence of both partition functions and the stimulated emission factor into the 
single factor (To/T)™. The exact temperature dependence of Boltzmann factors is retained. If the 
parameter m is available, use of (376) offers decreased computational overhead relative to (374). 

Pressure-broadened halfwidths o are also provided at the reference temperature and pressure, 
apfl = ap(po,To). hitran literature refers to broadening features are as air-broadened (rather 
than pressure-broadened) because all parameters in the hitran database are normalized to Earth's 
natural isotopic mixing ratios. Thus while the broadening process is collision-broadening, the 
tabulated line width refers to collision-broadening by air of Earth's isotopic composition, hence 
the label, hitran actually provides two halfwidths, the self -broadened halfwidth ag.o and the air- 
broadened or foreign-broadened halfwidth af^. Air-broadening accounts for collision-broadening 
by foreign molecules, i.e., all molecules except the species undergoing radiative transition. Thus 
air-broadening is, to a good approximation, due nitrogen and oxygen molecules. The other form of 
broadening, self-broadening refers to collision-broadening by molecules of the species undergoing 
the radiative transition. Because reasonances may develop between members of identical species, 
in general as,o 7^ ^f.o- The pressure-broadened halfwidth of transitions of radiatively active species 
A at arbitrary pressure and temperature scales as 

(^pip,T) = (^) (- — ^af.o + — as,o) (377) 
V^y V Po Po / 

where n is an empirical fitting parameter supplied by hitran and pA is the partial pressure of 
species A. Clearly self -broadening effects are appreciable only for species with signifcant partial 
pressures, i.e., O2 and N2. For minor trace gases, pa ^ P and pA <^ po so p — pa ~ p and 
Pa/Po ~ 0. For such gases, (377) simplifies to 

«(p,T) ^ ao-(^] (378) 
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which has no dependence on self -broadening. 

Collision-broadening may also cause a pressure-shift of the line transition frequency away from 
the tabulated line center frequency z>o = t'(po) to a shifted frequency D* = HITRAN supplies 
a parameter 5(po) in COS wavenumbers per atmosphere (cm^^ atm^^) with which to calculate the 
hne center frequency change due to pressure-shifting. 

z>*(p) = z>o + 5(po)- (379) 
Po 

The prescriptions for adjusting line centers (379) and half-widths (377) from (po, To) to to 
arbitrary p and T should be used to determine the corrected line shape profile of each transition 
considered. In the lower atmosphere, application of (379) and (377) leads to a corrected Lorentzian 
profile (280) 

1 «l(p,T) 



n {z/ - [i^o +p5(Po)/Po]}^ + aL(p,T)2 



<^>i^{u,uo,p,T) = -— — — (380) 

71 [u — u*)'' + ai^{p,T)^ 

Of course this section has only discussed application of hitran database parameters to atmo- 
spheres in LTE conditions. Application of HITRAN to nonlocal thermodynamic equilibrium condi- 
tions is discussed in Gamache and Rothman [1992]. 

It is necessary to compute statistics of 5*^ for use in narrow band models (§6). There is no reason 
not to apply the full approach of (374) when computing these statistics. Certain well-documented 
narrow band models use more approximate forms appropriate for specific applications. Briegleb 
[1992] uses 

S,{T) = Si,o(To) J g_hci>,/kTo 1 _ e-hc^.,/kTo (3^1) 
which is a hybrid of (374) and (376) with m = 1.5. 



5.6 Line-By-Line IModels 
5.6.1 Literature 

The classic paper on the water vapor continuum is Clough et al. [1989]. Fu and Liou [1992] com- 
pare a correlated-A; method to line-by-line model results. Ellingson and Wiscombe [1996] describe 
the a field experiment to directly compared measured and modeled spectral radiances to both band 
and line-by-line models. Crisp [1997] uses a line-by-line model to determine atmospheric solar 
absorption. Sparks [1997] presents an economical algorithm for selecting variable wavelength 
grid resolution such that absorption coefficients may be computed to a given level of accuracy. 
Mlawer et al. [1997] compare a correlated-A; model (rrtm) to the well-known line-by-line model 
LBLRTM. Moncet and Clough [1997] describe a fully scattering line-by-line model. Walden et al. 
[1997] presents test-case spectra for evaluating line-by-line radiative transfer models in cold and 
dry atmospheres. Vogelmann et al. [1998] use observations to constrain magnitude of any non- 
Lorentzian continuum in the near-infrared. Partain et al. [2000] discuss the equivalance theorem. 
Quine and Drummond [2001] present an algorithm that computes absorption coefficients to a spec- 
ified error tolerance by using a pre-computed lookup table of where interpolation is appropriate. 
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6 Band Models 

Band models, also called narrow band models, discretize the radiative transfer equation intervals 
for which the line statistics and the Planck function are relatively constant. The literature describing 
these models extends back to Goody (1957). Band models have traditionally been applied to 
thermal source functions (187a)-(187b), but they may also be applied to the solar spectral region 
[Zender et ah, 1997] (M0DTRAN3) if additional assumptions are made. 

6.1 Generic 

The following presentation assumes that the absorption path is homogeneous, i.e., at constant tem- 
perature and pressure. Important corrections to these assumptions are necessary in inhomogeneous 
atmospheres. These corrections are discussed in §6.3.2. 

This discussion makes use of an arbitrary frequency interval Az/ which represents the dis- 
cretization interval of narrow band approximation. It is important to remember that Az/ is best 
determined empirically by comparison of the narrow band approximation to line-by-line approx- 
imations. Typically, Az/ is 5-10 cm^^. Kiehl and Ramanathan [1983] showed Az/ = 5cm~^ 
bands are optimal for CO2. Briegleb [1992] uses Az/ = 5cm"^ for CO2, O3, CH4, N2O, but 
Az/ = lOcm-i for H2O. Kiehl [1997] recommends Az/ = Scm"^ for CO2, lOcm-^ for H2O, 
5-10 cm"^ for O3, and 5 cm"^ for all other trace gases. 

6.1.1 Beam Transmittance 

The gaseous absorption optical depth ra(z/) is the product of the spectrally resolved molecular 
cross-section a{v) (535a) and the absorber path U (269) 

ra(z/) = a{u)U (382) 

The transmittance T(z/) between two points in a homogeneous atmosphere is the negative expo- 
nential of Ta(z/) (182), 

T(z/) = e-"-^'^) 

= e-"(")^ (383) 

Note that T(z/) is the spectrally and directionally resolved transmittance because (a) it pertains to a 
monochromatic frequency interval and (b) it applies to radiances, not to irradiances (which require 
an additional angular integration). Narrow band models are based on the mean transmittance T^^, 
of a narrow but finite frequency interval Az/ between z/q — Az//2 and z/q + Az//2. 

Az/ J^^ 

= / e'^^dz/ (384) 



where z/q is the central frequency of the interval in question 
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Using (271), we rewrite (384) as 

Tau 



Au 



exp[—a{h')U] du 



If the lines are Lorentzian (i.e., pressure-broadened) then the transmittance may be written in terms 
of the line strength and absorber amount using (273a) and (305) 



a{u)U 



SaiJJ 



Thus (384) becomes 



Tau = ^ I exp 



Sai^U 



(385) 



(386) 



6.1.2 Beam Absorptance 

In analogy to T (181) we define the monochromatic beam absorptance A 

A = 1-T 

= 1 — exp[— a(z/)f/] 



(387) 



where a{u) is the absorption cross-section and U is the absorber path. As described in §6.1.1, we 
refer to A as the beam absorptance becuase it pertains to a monochromatic frequency, and refers 
to absorption of radiance, not irradiance. 

The band absorptance Aau is the complement of the band transmittance 



Aau 



1-Tau 
1 

1 



Adu 



Au 



1 - e di^ 



Au 



- — / 1 — exp[— a(z/)f/l du 

^^Jau 



(388) 



The relation between the band absorptance and band transmittance is the same as the relation be- 
tween the monochromatic beam absorptance and the monochromatic beam transmittance (387). 
Aau measures the mean absorptance within a finite frequency (or wavelength) range comprising 
many lines. Aau does not include any contribution from outside the Az/ range. Thus Aau may 
neglect contribution from the far wings of some lines in the band. Moreover, it is important to re- 
member that energy absorption follows the band absorptance only if the energy distribution within 
the band is close to linear. 

If the line shape is Lorentzian (280), then (386) applies so that 



Aau 



1 
A^ 



1 — exp 



Au 



Sa-JJ 



dv 



(389) 
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6.1.3 Equivalent Width 

The spectrally integrated monochromatic beam absorptance of a line is called its equivalent width 

/ + 00 
Adv 
-oo 

/ + 00 
1 - exp[-a{u)U] dz/ (390) 
-oo 

Thus the absorber amount must be specified in order to determine the equivalent width. The lower 
hmit of integration, — oo, is actually an approximation which is mathematically convenient to retain 
so that (390) is analytically integrable for the important hne shape functions. Recall that we are 
working with a frequency coordinate v which is defined to be z/ = at line center (306). In 
practice the error caused by using the analytic expressions resulting from rather than those 
resulting from Jl^^ is negligible except for lines in the microwave. 

The name "equivalent width" reminds us that W is the width of a completely saturated rectan- 
gular line profile that has the same total absorptance. The relation between U and W{U) is called 
the curve of growth. The curve of growth of a gas can be measured in the laboratory. The inter- 
pretation of the curve of growth played a very important role in advancing our understanding and 
representation of gaseous absorption. 

If the line shape is Lorentzian (280), then (389) applies 

r+oo 

W = 1 - exp 



Sai,U 



(391) 



6.1.4 Mean Absorptance 

Although the spectrally resolved absorptance (387) and the band absorptance (388) are dimension- 
less, the equivalent width (390) has units of frequency (or wavelength). It is convenient to define a 
dimensionless mean absorptance for a transition line A 



1 



+00 



A = ^ I Adv 



1 
1 
W 

T 



-oo 

+00 



oo 



1 — exp[— «(z/)[/] dz/ 

(392) 



where 5 is the mean line spacing between adjacent lines in a given band. Note that A/^y (388) 
differs from A (392). A accounts for the absorptance due to a single line over the entire spectrum 
and renormalizes that to a mean absorptance within a specified frequency (or wavelength) range 
(5). In contrast to A/\j,, A does not neglect any absorption in the far wings of a line. Therefore 
the sum of the average absorptances of all the lines centered within Az/ may exceed, by a small 
amount, the band absorptance computed from (388) 



J]A>^a. (393) 



i=l 
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This difference should make it clear that A is much more closely related to W than to Aai^ 
For lines with Lorentzian profiles (280), W and A are 



W = A6 



exp 



Sai^U 



Tl{i"' + ai) 

We shall rewrite (394) in terms of three dimensionless variables, x, y, and u 



X = 


u/S 


V = 


6x 


dx = 




du = 


6 dx 


y = 




CtL = 


5y 




SU 


u = 








SU 




n 





This change of variables maps u G [— oo, +oo] to x G [— oo, +oo]. We obtain 



A 



oo 
oo 



— oo 
oo 



oo 

+ 00 



— exp 
i — exp 
1 — exp 
1 — exp 



-{2a-Lu) 



2a'lu 



52(a;2 + y^) 
_2fyhi_ 
'P{x^ + y^) 

,2 1 



S^x'^ + b'^y'^ 
dx 

dx 



{bdx) 



2uy' 



x^ + y^ 



dx 



(394) 



(395a) 



(395b) 
(395c) 



(396) 



The solution to this definite integral may be written in terms of modified Bessel functions of the 
first kind ly{u) (see §10.4) 

A = 2Tiyue-''[Io{u) + h{u)] 

= 2'KyL{u) 

= 2nai^5-^L{u) 
W = 2naLL{u) (397) 

where L{u) is the Ladenburg and Reiche function [Ladenburg and Reiche, 1913]. The dimension- 
less optical path u (395) is therefore a key parameter in determining the gaseous absorptance. Two 
important limiting cases of (397) are u <^ 1 and 1. 

For small optical paths, SU <^ 1 (269), or, equivalently, u <^ 1 (395c). This is the weak-line 
limit. In this limit the exponential in (394) or (396) may be replaced by the first two terms in its 
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Taylor series expansion. Thus all line shapes have the same weak-line limit 



/ + 00 
a{u)Udu 
-oo 



oo 

/+00 
-oo 

^ SU (398) 

where we have use the generic normalization property of the line shape profile (283) in the final 
step. 

For large optical paths, SU ^ 1 (269), or, equivalently, m ^ 1 (395c). This is the strong line 
limit. To examine this limit we first note that the line HWHM is much smaller than the frequencies 
where line absorptance is strong, i.e., ai^ <^ v'm (394). Equivalently, y <^x (395a)-(395b) so that 
y"^ may be neglected relative to in the denominator of (396) 

/ + 00 
1 — exp(— 2ny^/x^) dx 
-oo 

One further simplification is possible. Both terms in the integrand are symmetric about the origin 
so we may consider only positive x if we double the value of the integral. 

A = 2 1 - exp{-2uy^/x^) dx (399) 
Jo 

The change of variables z = 2uy'^x~'^ maps x G (0, +oo) to z E {+oo, 0) 

z = 2uy'^x-^ (400a) 
X = {2ufl'^yz^l'^ (400b) 
dx = {2ufl''y{-\)z~^^'' dz 

= -2-^^\^/^yz~'^^^ dz (400c) 

Substituting this into (399) leads to 

A ^ 2 r (l-e-^)(-2-i/V/2y)^-3/Mz 

^+oo 

r" + 00 



2'/^u'/'y / {l-e-')z-^/'dz 
Jo 

/+00 
^-3/2 _ ^-3/2g-. 



The first term in the integrand of (401) is directly integrable J z"^^"^ dz = —2z~^l'^. The second 
term in the integrand of (401) is the complete gamma function r(— 1/2) =. Section ?? describes 
the properties of gamma functions'^. 



'^This section is currently in http : / /dust . ess . uci . edu/f acts/aer/aer .pdf . 
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Combining these results we find that the mean absorptance of an isolated Lorentz line in the 
strong-line limit (401) reduces to 

A ^ vV2u i -Iz-^l^ 



- (-2v/S) 



= ?/V2ti(oo + 2A/7r) 

= 2yV2jtu (401) 
fxm: The first term must vanish, but how? Substituting (395b)-(395c) into (401) 



S J V Sttol 

W = A6 = 2^/SUa[ (402) 

In the strong line limit we see that the mean absorptance A of an isolated line increases only 
as the square -root of the mass path. Physically, the strong line limit is approached as the line 
core becomes saturated and any additional absorption must occur in the line wings. Put another 
way, transition lines do not obey the exponential extinction law on which our solutions to the 
radiative transfer equation are based. One important consequence of this result is that complicated 
gaseous spectra must either be decomposed into a multitude of monochromatic intervals, each 
narrow enough to resolve a small portion of a transition line, or some new statistical means must 
be developed which correctly represents line absorption in both the weak line and strong line limits. 

In summary, the equivalent width W (390) of an isolated spectral line behaves distinctly dif- 
ferently in the two limits (398) and (402) 

7r f SU : Weak-line limit ^^^^^ 

^ = ^^ =\2VSU^^ : Strong-line limit ^^^^^ 



6.2 Line Distributions 

Inspection of realistic gaseous absorption spectra reveals that line strengths in complex bands oc- 
cupy a seemingly continuous distribution space, with line strengths S varying over many orders of 
magnitude in a single band. A key point discussed further in Goody and Yung [1989] is that the 
variabilities in line spacing and in line widths within a band are negligible (and usually order of 
magnitude smaller) in comparison to the dynamic range of S. A band containing a suitably large 
number of lines, therefore, may be amenable to the approximation that the line strength distribu- 
tion may be represented by a continuous function of S. The line strength distribution function 
p{S) is the probability that a line in a given spectral region will have a line strength between S 
and S + dS. The function p{S) must be correctly normalized so that probability of a line having a 
finite, positive strength is unity 

/•oo 

/ p{S)dS = 1 (404) 
Jo 

A number of functional forms for ^(5*) have been proposed. 
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6.2.1 Line Strength Distributions 

Goody [1952] proposed the exponential line strength distribution, now also known as the Goody 
distribution 

p{S) = S-^e-'^/^ (405) 

where S" is a constant which, in the next section, we show to be equal to the mean line intensity. 
A prime advantage of (405) is its arithmetic tractability. The zeroth and first moments of (405) are 
both solvable analytically. The following sections use this property to demonstrate the absorptive 
characteristics of line distributions. However, more complex distributions can improve upon the 
exponential distribution (405) by better representing the observed line shape distribution of many 
important species in Earth's atmosphere, such as H2O. 

Malkmus [1967] proposed an improved, albeit more complex, analytic form for the line strength 
distribution function. First we present a simplified, approximate version of this PDF which illus- 
trates the essential differences between the Malkmus distribution and simpler line strength distri- 
butions. 

p{S) = S-^e-^/^ (406) 

The difference with (405) is the prefactor has changed from S^^ to S^^. The S^^ dependence 
increases the number of weaker lines relative to stronger lines, which is an improvement over (405) 
and simpler distributions (such as the Elsasser or Godson distributions). The Malkmus distribution 
(406) is now perhaps the most commonly used line distribution function. 

The astute reader will note that (406) is not normalizable on the interval [0, +00) (cf. §6.2.2). A 
key point of the Malkmus distribution, therefore, is the truncation or tapering of p{S) so that statis- 
tics of the resulting p(S) are well-behaved. A straightforward line strength distribution function 
with most of the desired properties is the "truncated" distribution which is non-zero only between 
a fixed maximum and minimum line strength, Sm and Smin, respectively. By convention, Smm is 
deprecated in favor of the ratio Rs between the maximum and minimum line strengths considered 
in the band 

5'mm = Sm/Rs 

Rs = Su/Smin (407) 

The line strength ratio Rs measures the dynamic range of line strengths considered in a band, 
which can be quite large, e.g., greater than 10^ for H2O bands. With this convention, 

S < Sm/Rs 

Sm/Rs <S<Sm (408) 
S > Sm 

Godson [1953] appears to have been the first to examine (408), although it is closely related to 
the full Malkmus distribution. The normalization of (408) is demonstrated in §6.2.2. Clearly 
(408) contains no contributions from lines weaker than Sm/Rs or stronger than Sm- The main 
disadvantage to (408) is that it is discontinuous, and thus more difficult to treat computationally. 



p{S) 



(S In Rs 



u 

-1 
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Malkmus [1967] had the insight to identify the following continuous function 

= sh^lS 

We shall call (409) the full Malkmus line strength distribution. In practice, (409) is only applied 
in the limit Rs — )■ oo. Rather than simply truncating (406) so that p(S) is non-zero only between 
some lower and upper bound (408), Malkmus [1967] developed (409) because it has many useful 
properties: 

1. p{S) is continuous 

2. p(S) has S^^ dependence over the bandwidth of interest 

3. p{S) approaches (408) 

4. p(S) is analytically integrable 

Following Malkmus [1967], let us ellucidate the behavior of (409) in important limits. For 
large Rs and Sm/Rs <ti S <ti Sm the bracketed term in (409) approaches unity so that p{S) k, 
{SlnRsY^, i.e., the behavior is identical to (408). For S > Sm, p{S) ^ {S\nRs)-^e-^/^''' . 
Thus for S ^ Sm, piS) is non-zero but much less than what extrapolating {S In Rs)~^ into this 
region would yield. Finally, for large Sm and S <^ Sm/Rs, piS) ^ {Rs)~^{Rs — ^)/Sm- To 
summarize, the limits of (409) are as follows 



{\nRs)-\Rs-l)/SM 
PiS) = { (SlnRs)-' 
{SlnRs)-^e-^/^'' 



S ^ Sm/Rs 

Sm/Rs < ^ < (410) 
S > 5a/ 



Thus (409) behaves like (408) everywhere except near the truncation points Sm/Rs and Sm- Be- 
yond these points, (409) is sharply tapered so that lines in these regions do not influence the statis- 
tics of the line strength distribution very much. 

6.2.2 Normalization 

The normalization properties of the line strength distribution functions will now be shown. The 
exponential distribution (405) is easily shown to be normalized 

S-^e-^/^dS = S-^\-Se-^/^V 

I J 

= -0-(-l) = l (411) 

For heuristic purposes we first demonstrate that the approximate Malkmus distribution (406) is 
not normalizeable despite its apparent simplicity. With the substitution x = S/S we have 

POO 

S~^e-^/^dS = / x-^e-'^dx (412) 



The RHS resembles the complete gamma function of zero, r(0) (??). However, r(0) is indeter- 
minate, a singularity between negative and positive infinities much like tan ^. Thus (406) is not 
normalizeable. 
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Next we consider the truncated Malkmus distribution (410). We integrate over the non-zero 
region, Sm/Rs to Sm, to demonstrate its normalization properties 

"Sm 



/ {SlnRsr'dS = (Ini?5)-Mln5]|;/^^ 

JSm/Rs 



= {lnRs)-'[lnSM-HSM/Rs)] 

= i\nRs)-'[\nSM-\nSM + \nRs] 

= 1 (413) 

It is more challenging to prove that the full Malkmus distribution (409) is normalized. To begin, 
we simplify p{S) (409) with the substituions a = 5^/, b = RsS^^, and thus Rs = b/a. 



POO POO 

/ p{S)dS = (InRs)-^ S-^e-"'"^ -e-^^)dS 
Jo Jo 



(414) 



We must show, therefore, that this rather complex definite integral always equals ln(6/a) = In Rs- 
Since both terms in the integrand are of the same form, we shall explicitly evaluate S^^e^""^ before 
subtracting the analogous term involving b. These terms are amenable to integration by parts with 
the change of variables 



dv = S-^dS 
V = \nS (415) 



The result is 

r>00 



/»00 POO 

In Rs S^^e-^^dS = [e-"^ In S]"^ - In ^(-ae"''^) 
Jo Jo 

POO 

= [e-^^lnS]^ + a e-'^^ln5d5 
Jo 

The first term on the RHS contains a singularity since ln(0) is undefined. This term however, will 
disappear once the similar term arising from the integration of J S^^e~^^ dS is subtracted from 
(416). Now we change variables to isolate the parameter a from the dummy variable of integration. 
Letting x = aS, S = x/a, dx = a dS, and dS* = a^^ dx leads to 

POO POO 

In Rs ^-^e-'^^d^ = [e-''^\nS]'^ + a e"^ ln(x/a)a-Mx 
Jo Jo 

POO 

= [e""*^ In 5] ^ + / e'^'lnx - e'^'lnadx 
Jo 

POO 

= [e-''^\nS]'^ + / e-^'lnxda;-lna [-e-^]^ 
Jo 

POO 

= [e-''^\nS]^+ e-^'lna;da;-lna[-0- (-1)] 
Jo 

POO 

= [e-^^\nS]'^+ / e-^lnxdx-lna (416) 
Jo 
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There is no closed form solution to the second term on the RHS, / e~^lna;dx, which is very 
close in appearance to the original integral (414). However, this term will also disappear once 
the integration involving b is subtracted from (416). Multiplying both sides by (Ini?^)^^, and 
subtracting the integration involving b leads to 



+ / e In X dx — / e ^ In x dx — In a + In 6 

Jo Jo 
In 6 — In a ln(6/a) Ini?^ 

InRs ~ InRs ~ InRs 
= 1 (417) 

All the like terms in the numerator cancelled in the first step with the exception of In Rs- 



6.2.3 Mean Line Intensity 

The mean line intensity 5 of a line strength distribution p{S) is defined as 



oo 



S = / Sp{S)dS (418) 







5* is an important statistic of a line strength distribution and any realistic line strength distribution 
should adequately predict S when compared to observed line strengths. Thus many line strength 
distributions, e.g., the exponential distribution (405) contain S in their definition. The claim that S 
is, in fact, the mean strength of these distributions will now be proved. 
The mean line intensity of the exponential distribution (405) is 



/'OO POO Q 

/ Sp{S)dS = / te-'/'dS 
Jo Jo ^ 



1 

1 
S 



oo 



Se-'/^dS 







oo 



Sxe-^'Sdx 



oo 



S xe ^'dx 



5[-0-0-(0-l)] 

S (419) 



This demonstrates S is the mean line strength of the exponential distribution. 
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The mean line intensity of the approximate Malkmus distribution, (406) is 

S 



Sp{S) dS 







S 



dS 



e"^/^ dS 

, - oo 





-oa/S _ 



-0/5N 



(420) 



This proves that 5* is also the mean line strength in the approximate Malkmus line distribution 
function (406). 

The mean line intensity of the full Malkmus distribution is derived directly from (409) 



Sp{S) dS 



S X 



-RsS/Sm 



SlnR. 



dS 



,-1 



(In i?5 
1 

Sm 
Sm , 
Rs-i 

RslnRs 



^-S/Sm _ ^~RsS/Sm 



Sm_ r^S/Sm 



Rs 

[-0 + -(-!)- (^5' 
Sm Rs - 



1 - 



S 



M 



(421) 



This expression for S does not converge in the hmit as i?5 — )■ 00. Thus the mean line strength of 
the Malkmus distribution (409) depends on the parameter Rs (??), the ratio between the strongest 
and weakest lines in a given band. It is clear from (421) that S decreases as Rs increases. Ade- 
quately representing the importance of weak lines is one of the chief advantages of the Malkmus 
distribution relative to other distributions. 

The maximum line strength intensity Sm may be expressed in terms of the mean line strength 
S by inverting (421) 



S. 



R^\nR. 



M 



R, 



1 



S 



{All) 



6.2.4 Mean Absorptance of Line Distribution 

This section develops the concept of mean absorptance of a band of non- overlapping lines, whereas 
§6.1.4 derived the mean absorptance A of an isolated Lorentz line. Since the lines are assumed to 
be non-overlapping, we may weight the contribution of each line to the mean absorptance (392) by 
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W = A6 



p{S) 



+00 



Adu dS 



oo 

+ 00 



poo p-i-oc 

/ p{S) / 1 -exp[-S$f/]dz/d^ 

Jo J~oo 



(423) 



Inserting the exponential line strength distribution function (405) for p{S) in (423) and inter- 
changing orders of integration yields 



W = A6 



+ 00 POO 



+00 



e - exp 



1 -exp[-5$f/] dudS 
S 



-oo ^0 

+ 00 POO 



s 



S<^U dSdu 



-s/s 



-oo Jo 

+ 00 



exp (-5(^-1 + f/$)) d^d 



oo 

+ 00 



5e-^/^ + {S^^ + f/$)-i exp {-S{S~^ + f/$)) 

/ + 00 
[0 + - (-5 X 1 + {S-^ + t/$)-i X 1)] du 
-oo 

/ + 00 
S ~ {S-^ + U<^)-^ du 
-oo 



du 



+ 00 



oo 

+ 00 



s 



s 



du 



1 



1 



1 + su<^ 
1 + 5f/$ - 1 



— oo 

+ 00 



1 + SU<^ 



du 
du 



1 + SU<^ 



du 



(424) 



Thus the mean absorptance of the exponential line distribution depends straightforwardly on the 
line shape function. Substituting the Lorentz line shape $l(^) = + ^l)] (280) into (424) 

we obtain 



W = A6 



+ 00 



X 



SU 



du 



n 



1 



■oo 

+ 00 

■oo 



2 _L ^,2 n{u^+al)+SUa-[^ 



1 



dv 



dv 
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SUai^ 

SU ttL 



+00 



+00 



u"^ + al + n-'^SUaL 



(425) 



Letting a = a/c^l + T^'^SUai, = ai, a/1 + ^[//(TCaL) and noting that / (x^+a^) ^ dx 



a ^ tan ^ ( 



TT a 



+00 



+ 



dz/ 



tan 



+00 



5f/ 1 + 



5f/ 



SU 



w 

T 



A 



SU 



1 



SU 



-1/2 
-1/2 



1/2- 



L2 



(426) 



In addition to U, the mean band absorptance A depends on only two statistics of the band, namely 
5/5 and S/inaL). 

We recall that (426) is based on the idealized, continuous line strength distribution (405) rather 
than on real (observed) line strengths. Observed line strengths and positions are tabulated in the 
HITRAN database (§5.5.1). Therefore we may compute any desired statistics of the actual line 
strength distribution. Of particular interest is the observed mean absorptance of a band in both the 
weak and strong line hmits. Summing (398) and (402) over the N lines in the given band 



W = A6 



(427) 



The weak and strong line limits of the mean absorptance of the exponential line strength distribu- 
tion function (426) are 



lim W = AS^SU 

5C/<1 

lim W = A5^ SU 

5C/>1 



SU 



-1/2 



^TiSUaL 



(428a) 
(428b) 



We now require (428) to equal (427) in the appropriate limits. In the weak line limit 



N 



SU = N-^J2^iU 



i=l 



N 



iN6r'j2s. 



(429) 



i=l 
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In the strong line limit 



N 



i=l 

AU 



N \ 

1=1 J 



1 N' 



2 



N 



1 /V2 / ^ 

= S^X^^v^CL,, 



, i=l 
2 



S 



2 



\i=l ' i=l / 



The tractable analytic form of (426) led to its adoption as the dominant formulation for band 
models for many years. The mean band absorptances A of other line strength distributions have 
forms similar to (426). 

Inserting the Malkmus line strength distribution function (409) for p{S) in (423) and inter- 
changing orders of integration yields 

/•oo p-S/Sm _ p-RsS/Sm 
W = A5 = / — — / 1 -exp(-5$f/)dz/d5 

Jo ^ 111 J-oo 

"+00 POO 



In 



1 p-rco poo 1 

L_ / / i-(e-^/^M - e-^^^/^M)(l - e-^*^) dS di 

J-oo Jo ^ 



^ I-+00 i^oo ^-S/Sm q-SRs/Sm Q-Sm+Rs/SM) q-S{W+1/Sm) 

H t; t; dS du 



In Rs J —oo J ^ ^ ^ 



The first two terms and the last two terms in the integrand are both of the form S~^{e~°'^ — e"'''^). 
In proving the Malkmus distribution is correctly normalized (414)-(417) we showed the value of 
this type of definite integral is ln(6/a). For the first two terms, a = Sj.^ and b = Rs/Sm- For 
the final two terms, a = $[/ + Rs/Sm and b = + 1/Sm- Therefore the integration over S is 
complete and results in 



1 , ..f ^u+i/s. 



W = A6 = — I ln/?5 + ln( ' 1 du 

In Rs J. ^ \m + Rs/Sm ^ 

1 /■+°° , „ , M + Sm^U 



InRsJ.^ \Rs + Sm'^U 



1 



lni?5 



1 \ I + Sm'^U 



+00 

In Rs + In 

■oo 



RsJ 1 + Sm^U/Rs 



du 
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-— / lni?s + ln — — — — — l-lni?sdz/ 
In ^5 i-oo V 1 + Sm^U /Rs , 



1 ''"iJ'l'^fL Id. (431) 



"+00 

The mean absorptance of the Malkmus distribution also depends simply on the line shape $(z/). 
The mean absorptance of the exponential line distribution (424) is analogous to the Malkmus result 
(431). Note that the former is a rational function of the mean line intensity S, while the latter de- 
pends logarithmically on S. Moreover, the mean absorptance of the Malkmus distribution depends 
on the additional parameter Rs, the ratio between the maximum and minimum line strengths with 
the band. The dependence is problematic in that (431) is not self-evidently convergent as Rs — )■ oo. 
We now show that for the most important line shapes, ^(431) does converge in this limit. 
Substituting the Lorentz line shape = oii,/[n{i''^ + a^] (280) into (424) we obtain 

W = AS = ^ I In I 1 du 




SmUulRs' + T^{y^ + al) 

u"^ + al + SMUaj^/n 
u^ + al + SMUai^/{nRs) 

z/ + a£ H ] — In [ u + ai-\ — — ) du (432) 



du 



du 



We simplify (432) by making the substitutions = oI+SmU aL/Jt and cP = o^+SmU ai,/{nRs 
or, equivalenty 



aLA/l + ^ (433a) 



d = a^dl + ^^ (433b) 



With these definitions, (432) becomes 



1 



+00 



W = A6 = —— / ln(z/' + c') - ln(z/' + d') dz/ (434) 

This is the difference of two identical, standard integrals [e.g., Gradshteyn and Ryzhik, 1965, 
p. 205] and reduces to elementary functions 



/ ln(z/^ + a^) du = u ln(z/^ + a^) - 2z/ + 2z/ tan"^ - 
J a 



(435) 



It is straightforward to differentiate (435) and verify the result. Using (435) to complete the fre- 
quency integration in (434) results in 

_ r 1/ +00 

W = A6 = (InRs)-^ z/ln(z/2 + c2)-2z/ + 2ctan-^--z/ln(z/2 + rf2) + 2z/-2dtan-^ 



c d. 



(In/?, 



s) 



2c tan 2a tan — 

c d- 



+00 
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The mutual cancellation of the 2z/ terms is obvious, but the equivalence of the logarithmic terms 
may not be self-evident. Consider, however, that both c and d (433) are finite-valued as — )■ ±00, 
and it is clear that these terms must cancel. Evaluating the remaining terms leads to 



W = A6 = {In Rs)'^ {2c^~2d'^ 



2 2 
(In Rsy^{nc -nd + Kc- nd) 

(In Rs)~^2n{c-d) 



Although (436) does convey the behavior of the mean absorptance of the Malkmus band model 
with Lorentzian lines, the presence of the Sm term and the Rs term is undesirable because their 
values are interrelated and depend on the specific gaseous absorption band being considered. It 
would be more convenient to express the mean absorptance in terms of universal statistics appli- 
cable to all band models, so that A and W may be easily intercompared among band models. 
The mean line strength S (418) appears in both the exponential (405) and Godson line strength 
distributions. We use (422) to replace Sm by S in (436) 



In Rs \ y ttol -Rs - 1 V ^"l-Rs Rs - ^ 



The weak fine hmit of the mean absorptance of the Malkmus line strength distribution function 
(437) is obtained by expanding the radical into a Taylor series and retaining only the first two terms 

lima-^i a/1 + X ^ 1 + x/2 

- 2Tia^ f SU RslnRs , SU InR. 
lim W = Ao ^ 1 H 1 



sU't:i In-Rs \ 27raL Rs — 1 2'Ka-L Rs — 1 

27raL SU\nRs(Rs-l) 

— X 



\nRs 2nai^{Rs~l) 
SU (438) 



which is identical to (428a). The strong line limit of (437) is 



lim W = A5 

su->i 




2 



rTlai^SU\ 
\h[Ri) 
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= 2^/%ai^SUx , (439) 
^{Rs-l)\nRs 

Let us encapsulate the i?5 -dependence into the function x{Rs) and and also present x^{Rs) for 
future use. 

x{Rs) = , ^ = (440a) 

v/(i?5 -1) In i?5 



2 



^-l)(v/i!^-l) 

(i?5-l)lni?5 
(/R^-l)(v/i!^-l) 

^-l)(v^+l)ln/?5 



where we have made use of the factoring Rs — 1 = (y/Rs — ^){VRs + !)• Rewriting (439) with 
x{Rs) (440a) 

lim W = A6 ^ 2^/naj^SUx{Rs) (441) 

5C/>1 

The strong line limit of the Malkmus distribution of Lorentzian lines (441), and the strong line 
limit of an individual Lorentzian line (403), differ by the factor y/nxiRs)- 

The advantage of encapsulating the i?5 -dependence in the xi^s) function is that the weak and 
strong line limits of the Malkmus distribution may now be cast in exactly the same notation as the 
exponential distribution (428). First introduce a transformed mean line intensity S and mean line 
spacing 5 which are defined in terms of the line strength ratio of the Malkmus distribution 



S 

nx"" 



S = S = nSx^ (442a) 



I = — ^ 5 = TT^x^ (442b) 

The definitions (442) ensure that S/5 = S/5. In terms of the original mean line strength and 
spacing, and the transformed line strength and line spacing, the weak and strong line limits of the 
mean absorptance A of the Malkmus line strength distribution apphed to Lorentzian hnes (437) is 

lim W = A5^ ^ = — (443a) 

sr/<i 5 § 



lim W = A6^ 2^'!iai^SUx{Rs) = '^y a^SU (443b) 

S'C/>1 

We now rewrite the mean absorptance (437) in terms of the transformed parameters S and 6 



6 BAND MODELS 



113 



(442), making use of (440b) 



W 



A 



6 In Rs 



1 2naL 




nx^SU \nRs 



^-1 



6\nR, 




Rs-l)iVR~s + l) 



^-1 



1 H X — =^ X 

ctL (/R^ +1) In i?5 



/i?^-i)(y:R^+i) 




The asymmetric nature of the optical path terms under the radicals in (444) is clear. For large 
Rs, the first absorptance term approaches a finite-valued function of S, U, and q;l. The second 
absorptance term, on the other hand, approaches unity. Thus the mean absorptance A asymptotes 
to a finite value as -R5 — )■ 00. Formally, this amounts to taking the limit as ^ 1 



W 
T 



A 




VT+0 



(444) 



To obtain W and A in terms of S and 5 we use (440a) in (444) 



W 



A 



2aL X 




(445) 



The mean absorptance of the Malkmus-Lorentz distribution A (445) is thus a straightforward func- 
tion of S, U, a-L, and the parameter x (440a). Equation (445) is valid only in the limit y/Rs 3> 1. 
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It is more computationally economical to choose a single value of x than to determine a new 
value for each band considered. The traditional approach to this has been to choose x such that the 
mean absorptance of the Malkmus distribution (445) equals the mean absorptance of the exponen- 
tial distribution (426) in the weak and strong line limits. The weak and strong line limits of (445) 
are 



Jim W = A5^ 27raLX^ 1 + 



SU 



1 



SU 



(446a) 



Jim W = A6 ^ 2'Kai,x 

5i7>l 




2xV'!^SUaL 



(446b) 



Equations (446a) and (428a) are an identity and yield no information on the value of x{Rs)- Equat- 
ing (446b) to (428b) yields 



2x = 1 
X = 1/2 



(447) 



Thus the value X = \ causes (446) to agree with (428). Substituting = 1/4 into (445) yields the 
traditional form of the mean absorptance of the Malkmus line strength distribution for Lorentzian 
lines 



W_ 
5 



A 




(448) 



As mentioned above, the choice x^ = \ (447) implies a specific choice of Rs- Starting from 
(440b) 



1 



ms + l)\nRs 
InRs 



1 
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(449) 



which is a transcendental equation in i?^. The solution of (449) appears to be imaginary, which is 
inconsistent with the assumption that ^ 1 (444). 
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6.2.5 Transmittance 



The band transmittance T^u of the exponential distribution (405) is the negative exponential of 
(426) 
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Au 
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SU 
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The band transmittance T^u of the Malkmus distribution (409) is 



(450) 
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2(S/7raL)A. 
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(451) 



The band flux transmittance T^^^ is simplified by replacing the angular integration with a dif- 
fusivity factor, D. 
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(452) 



As discussed in §2.2.15, D = 1.66 gives the best agreement between (452) and models which 
integrate the RT equation directly 

Applying the Malkmus model to shortwave spectral regions is also possible if we redefine the 
diffusivity factor. For solar radiative transfer in clear skies most of the downwelling radiation 
follows the path of the direct solar beam. Thus the downwelling radiation traverses a mass path 
which increases with the cosine of the solar zenith angle 9. Defining yU = cos 9 we have 
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/i \7t:q;l/a 
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(453) 



The role of ^ in (453) is analogous to the role of D in (452). Note that both fi and D are approx- 
imations to the angular integral. Diffuse radiation, i.e., thermal emission, upwelling (scattered) 
solar radiation and downwelling solar radiation beneath optically thick scatterers, is better charac- 
terized by D = 1.66 than by /i. Direct downwelling solar radiation, on the other hand, is better 
characterized by /i than by D = 1.66. 



6.2.6 Multiplication Property 

It is observed that the spectrally resolved transmittance of a mixture of gases obeys the multipli- 
cation property. Let T(l, 2, ... A^) be the spectral transmittance of a mixture of A^ gases, and T{i) 
be the transmittance of the ith gas alone. 

N 

T{l,2,...N) = l[T{z) (454) 

i=l 

The usefulness of the multiplication property of transmittance cannot be overstated. It is basic to 
the construction of efficient wide band (emissivity) models. 
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Equation (454) may be integrated over a finite frequency interval. 

/ T(i,2,...iv)dz/ = / n^(^)^ 



N 

i=l 
N 

i=l 

N 



Tau{N) = [ l[T{t)di 

N 

= X{Tau{i) (455) 



i=l 



The final step is strictly valid only if the spectral transmittances of the N gases are uncorrected 
over the range of Az/. This will always be the case if the lines in these gases are randomly dis- 
tributed in the interval. Fortunately, this random band assumption is met by all important atmos- 
pheric bands. Henceforth we shall consider (455) valid unless otherwise specified. 



6.3 Transmission in Inhomogeneous Atmospheres 

Until now we have focused on methods designed to obtain accurate spectral mean absorptance and 
transmittance of one or many lines over homogeneous atmospheric paths, i.e., paths along which 
temperature, pressure, and gas mixing ratio do not vary. This is generally a reasonable assumpution 
only along horizontal atmospheric paths, i.e., at constant altitude. In realistic atmospheres, at least 
pressure and temperature will vary vertically, and cause corresponding changes in line strength and 
fine width (374). 



6.3.1 Constant mixing ratio 

Exact, analytic solutions for mean transmission through inhomogeneous atmospheres exist for 
Lorentzian lines in the idealized case where the gas mixing ratio is uniform, the atmosphere is 
isothermal at temperature T, and where, therefore, the half- width depends only on pressure p. Let 
q be the gas mass mixing ratio, and p be the atmospheric mass density so that the total gas mass 
concentration is qp. The Lorentzian line shape has a linear dependence on pressure (298). We shall 
assume hydrostatic equilibrium applies, 

dp = —pg dz 

dz = -{pg)-'dp (456) 

and use (456) to change variables to pressure coordinates. Referring to (271) and (273a), we 
may write the absorption optical depth through a uniformly mixed layer in hydrostatic equilibrium 
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between heights z and z' as 



S^j^{z")dU 
SM^")qp{z")dz" 
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'\2 



(457) 



where we have used the simple relationship U = qp/g to determine the mass path of a gas with 
constant mixing ratio q. 

Near the top of the atmosphere, the line shape approaches the Doppler shape $d (311). In an 
isothermal atmosphere ckd does not depend on height. The optical depth from level z to space is 



ra(^, z') 



Sq_ 
2%g 



In 



(458) 



6.3.2 van de Hulst- Curtis- Godson Approximation 

The general problem is to account for the effects of inhomogeneous mass, temperature, and pres- 
sure paths on transmission. We shall attempt to determine an scaled absorber path U, scaled 
temperature T, and scaled pressure p such that U, T, and p yield optimal results when employed 
in our machinery for homogeneous paths, e.g., (450), (451). We shall show that the scaled param- 
eters can be chosen to yield correct band transmittances in certain limits. 

The extra degree of freedom permitted by each parameter means that U, T, and p form a three- 
parameter scaling approximation. The researchers van de Hulst, Curtis, and Godson independently 
discovered the following scaling approximation, now known as the H-C-G approximation. First 
we note that absorption along a path in Earth's atmosphere is generally more affected by changes 
in p than by changes in T. This stands to reason since p changes by many orders of magnitude 
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from the surface to space, while T generally changes by less than 30%. Thus the usualy practice, 
which is often adeequate, is to set T equal to the mean temperature along a path, 

f = f (459) 

This leaves two parameters, U and p, to account for the effects of inhomogeneity on path absorp- 
tion. We first derive a single parameter scaling approximation based solely on U. This will prove 
useful afterwards in developing the two parameter H-C-G approximation. 

Consider an the absorption optical depth through an arbitrary inhomogeneous atmosphere 
whose mass absorption coefficient ^(p, T) (277) depends only on T and p. The optical path 
through this medium is (271) 

^{p,T)dU (460) 

A one-parameter scaling approximation is exact if ijj{p,T) can be separated into a frequency- 
dependent function ?](//) times a p-T dependent function T) [Goody and Yung, 1989, p. 224] 

V^(p,T) = $(p,T)r/(z/) 

$(p,T) = ij{p,T)/r]{u) (461) 
Applying this separation of variables (461) to (460) yields 

^$(j9,T)r/(z/) dU 



^(P.r)W.r) 1^ 
*(p,r) 'K") 

We see that ipip^T) is independent of the variable of integration, U , since both p and T are path- 
mean quantities. Therefore il){p^ T) may be placed outside the integrand 

= ^{p,T)U where 
if ^ f^iU (462) 

The scaled absorber amount, U, together with the absorption coefficient evaluated for the scaled 
temperature and pressure, exactly predicts the absorption optical path (460). 

The exact form of the scaling function T) is chosen to optimize the accuracy of (462). The 
definition (461) makes clear that T) is a weighting function that does not explicitly depend 
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on frequency, and yet must somehow account for the temperature dependence of all the lines with 
the band. Line strengths are highly temperature dependent (374). Considering first Doppler lines, 
we note that the Doppler line width q;d (310) has no pressure dependence, but does depend on 
^/T. Both strong and weak Doppler lines within a given band will therefore be well-represented 
by identifying the scaling function T) with the line strength S 

N 

$(p,T) = 5^Si(T) (463) 

i 

Thus T) has no explicit pressure dependence in this case. 

Chou and Arking [1980] present a more formalized one-parameter scaling based on the pre- 
sumption (appropriate for water vapor) that errors due to inhomogeneity are likely to bias line wing 
absorption more than line center absorption. Liou [1992], p. presents a detailed account of this 
scaling. 

The H-C-G approximation, sometimes called the C-G approximation, . . . 

6.4 Temperature Dependence 

SWNB evaluates the temperature dependence of the random band parameters. Let Tq be the tem- 
perature at which the line parameters (strengths, half- widths) are known. Typically Tq = 296 K, 
the reference temperature of the HITRAN database discussed in §5.5.1. The H-C-G approxima- 
tion accounts for the influence of temperature and pressure inhomogeneity along the path on the 
absorption (460). The scaling functions $(T) and \&(T) at an arbitrary temperature T are 

$(T) ^ 5^5,(T) /5^5,(To) (464) 

i ' i 

^(T) = A/'5.(T)aL(po,T)) / v/5,(To)aL(po,To)) (465) 



where S'j(T) is scaled from ^^(To) using (374) and ai^{j)Q,T) is scaled from aL(PO)To) using 
(377). 



The scaled absorber path U is 



The scaled pressure p is 



U = I q^{T)^ (466) 



p = i / q^{T)^^ (467) 

U J Po g 

Following Rodgers and Walshaw [1966], we fit the logarithms of $(T) and \E'(T) to quadratic 
polynomials in temperature 

In $(T) = a$(T - Tq) + 6$(T - Tq)^ (468a) 
In ^(T) = a^(T - Tq) + 6.t(T - To)^ (468b) 
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The linear term dominates (468a) since $(T) (463) depends most exponentially on T (374), i.e., 
^ b^. On the other hand the quadratic term is important in (468b) since \I/(T) (??) contains 
significant nonlinear behavior. The coefficients a$, 6$, a^, and are determined by least-squares 
fitting $ and \1' at discrete temperatures Tk over the temperature range of interest: 

^ _ - To) ln$(T,) - To)^ - E.(T, - Tp)^ In $(T,) - To)^ 

, _ Ekin - To? MjT,) Y.,{Tk - To)^ - - To) In $(T,) - Tp)^ 

E.(7^;^-To)2E.m-Tor-E.m-To)^E.(r.-To)3 ^ ^ 

and and 6^ are defined analogously. 

In practice SWNB adjusts all band parameters from the HiTRAN-standard temperature (296 K) 
to 250 K, which is near the mass-weighted mean temperature of Earth's atmosphere. This reduces 
any biases caused by non-quadratic temperature-dependent behavior of the logarithms of $(T) and 
-<^{T) (469). $(T) and -^\t) are evaluated exactly in the temperature range 180 < Tk < 320 K 
with a 1 K resolution and the results are fit using (469). The expressions (469) are generally 
accurate to within 1% over this 140 K range. 



6.5 Transmission in Spherical Atmospheres 

Until this point we have assumed plane-parellel boundary conditions. That is we have assumed the 
increase in path length due to off-vertical transmission is given by the secant of the solar zenith 
angle times the vertical path. Planetary atmospheres are, of course, spherical and the plane parallel 
approximation must be corrected for spherical effects when the zenith angle is large. Horizontal 
irradiances must be corrected if the zenith angle is greater than about 80°. 



6.5.1 Chapman Function 

The earliest quantitative computations of the absorber path to space at arbitrary angles through a 
spherical atmosphere were made by Sidney Chapman. Let H be the scale height of the species, so 
that 

n{z) = n(zo)e-("-"«)/^ (470) 

The scale height well-mixed gases in Earth's atmosphere is about 7 km. Let be the Earth's 
radius, and zq be the altitude at which the absorber path is to begin. With the assumption that the 
species is distributed exponentially with altitude (470), the number path N from height zp to space 
at a polar angle 9o reduces to 

pz=oo 

N{zo,9o) = / niz) 

J z=zo 

It can be shown that the number path (471) is a function of a dimensionless measure called the 
radius of curvature of the absorber distribution X 



1 - 



+ Zq) sin( 



-i/z 



(rffi + z) 



dz 



X = (re + zo)/if 



(471) 
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According to Huestis, X varies from about 300 to 1300 on Earth. The Chapman function Ch [e.g., 
Rees, 1989] is defined implicitly by the RHS of (471). It directly relates N{zo,9o) to the scale 
height H and local number concentration n(z) of a species via 

iV(zo,^^o) = Hn{zo)Ch{X,eo) (472) 
Note that Ch (472) depends only on X, not on the individual values of r^, zq, and H. 

7 Radiative Effects of Aerosols and Clouds 

Interaction of radiation and matter occurs on three scales. On the molecular scale, the scattering 
and absorption of radiation is described by the ambient concentration, scattering cross sections, 
reasonant frequencies, and quantum states of the molecule. Interaction of radiation with bulk 
phase matter (e.g., the sea surface) is described by the change in the index of refraction of the 
media (e.g., Snell's law) and macroscopic thermodynamic properties such as = aT"^ (117). In 
between these two scales is the interaction of radiation with atmospheric particulates, i.e., aerosols 
and hydrometeors. Particulates are, by definition, isolated bodies of condensed matter surrounded 
by a continuous media. We now turn our attention to describing the effects of particulates on the 
radiation field. 

7.1 Single Scattering Properties 

Three properties are required to exactly specify the radiative effects of particles. The properties 
define the total extinction (scattering plus absorption) due to the particle, the probability an inter- 
action results in absorption rather than scattering, and, finally, the angular distribution of scattered 
photons as a function of the incident angle. By convention, the properties defining the above at- 
tributes are usually specified as the extinction optical depth t^, the single scattering albedo w, and 
the asymmetry parameter g. These three parameters, Tg, w, and g, are known collectively as the 
single scattering properties of the particles. These properties depend in turn on the mass, size, and 
composition of the particle species. 

A particle's chemical composition determines its index of refraction}^ n. The index of refrac- 
tion is a dimensionless, complex number 

n = Uy. + irti (473) 

All the physics describing the fundamental electromagnetic properties of the material are consoli- 
dated into n. Uj- describes the scattering properties of the medium while rii describes the absorption 
properties of the medium. The literature disagrees about the sign convention for Ui. We choose to 
represent rii as positive-definite, i.e., ni > 0. This convention is employed by Bohren and Huffman 
[1983], among others. 

rij- and Ui are fundamental properties of matter and must be determined from laboratory studies 
[e.g., Hess et ah, 1998; Patterson et al, 1977; Patterson, 1981; Perovich and Govoni, 1991; Tang, 

'^Care must be taken not to confuse the complex index of refraction n (dimensionless) with the size distribution 
nil [#m^'^m^^]. The former is a function of wavelength for a given chemical composition, while the latter is a 
function of size for a given particle population. The meaning should be clear from the context. 
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Huffman, 1983, p. 12] 




V D 


= PF 


VB 
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V X H 


= Jf + 



1997; Volz, 1973; Bohren and Huffman, 1983], or from models constrained by field measurements 
[e.g., Colarco et ah, 2002; Sinyuk et ah, 2003]. The HITRAN database [Rothman et ah, 1998] 
contains a compilation of n for many aerosols of atmospheric interest. 

7.1.1 Maxwell Equations 

The macroscopic form of Maxwell's equations in Rationalized SI (MKSA) units is [Bohren and 



Coulomb's law (474a) 

(474b) 

Faradays's law (474c) 

— — Maxwell's generalization of Ampere's law (474d) 

dt 

where E is the electric field measured in Volts per meter [Vm^^], B is the magnetic induction 
measured in Tesla [T], and Jf is the/ree current density measured in Amperes per square meter 
[Am~^]. In this context, macroscopic refers to spatio-temporal averaging of the Maxwell equa- 
tions in a vacuum. Clearly (474) applies to microscopic spaces, since it is expressed in terms of 
differentials, though it does not apply at the smallest molecular scales. 
The electric displacement D and magnetic field H are defined by 

D = eoE + P (475a) 
H = — - M (475b) 

where P is the electric polarization (mean electric dipole moment per unit volume), M is the 
magnetization (mean magnetic dipole moment per unit volume). P is measured in [Cm~^], and 
M is measured in [A m~^]. The scalar symbols are the/ree charge density pp, measured in [C m^^], 
the permittivity of free space Eq, measured in [Fm^^], the permeability of free space po, measured 
in [Hm^]. 

The constitutive equations in SI units are [Bohren and Huffman, 1983, p. 13] 

Jf = ctE (476a) 
B = /iH (476b) 
P = £oxE (476c) 

Here a is the conductivity which is measured in Siemens [S], Amperes per volt [A V^^], or mhos 
per meter [mhox~^]. p is the permeability measured in [Hm^^]. x is the electric susceptibility 
measured in Farads per meter [Fm^^] or [C V^^ m^^]. Together, a, fi, and x known as the 
phenomenological coefficients — they define the medium under consideration. In a linear medium, 
these coefficients do not depend upon the electromagnetic field. In a homogeneous medium, these 
coefficients are independent of position. In an isotropic medium, these coefficients are independent 
of direction. The media we consider is assumed to be linear, homogeneous, and isotropic. 
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We shall now simplify the Maxwell equations (474) to the problem of scattering of radiation by 
matter described by the the constitutive equations (476). Subsituting (476c) into (475a) and (475a) 
into (474a), we obtain 

V-D = V- (£oE + £oxE) =£o(l + x)V-E = pF (477) 

Subsituting (476b) into (474b) leads to 

V ■ B = V ■ (/iH) =/iV-H = V- H = (478) 

Subsituting (476b) into (474c) leads to 

V X E + = (479) 

ot 

Subsituting (476c) into (475a), and then sustituting (476a) and (475a) into (474d), we obtain 

VxH = .E+^-^!^^^^^aE + e„(l + x)f (480) 
where e, the complex permittivity is defined by 

e = £o(l + x)+i- (481) 
(jj 

We assume the time-variation of the electro-magnetic field is harmonic, i.e., of the form E oc 
exp(— iwt) where uj is the angular frequency. The particular problem we wish to solve is the 
scattering of incident plane waves by a sphere. The harmonic plane wave definition (237) is 

Ec(r,t) = Eoexp[i(k • r - wt)] 

(482) 

where the subscript c denotes complex plane wave. With this assumption, the Maxwell equations 
for the medium (477)-(480) impose a dispersive relationship upon wave characteristics of (237). 

V ■ Ec = (483a) 

V ■ He = (483b) 

V X Ec = iw/iHc (483c) 

V X He = -iweEc (483d) 

Electrodynamics texts are often written in CGS (Gaussian) units. The Maxwell equations in all 
common units appear in Jackson [1975] p. 818. 

The coupled vector equations (483) will prove more tractable once they are reduced to coupled 
scalar equations. This will be accomplished by exercising (483) with some vector mathematics 
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(560). First we take the curl of equations (483c) and (483d): 

VxVxEc = Vx (iw/iHc 
= iwyuV X He 



iu;/i(— iweEc) 

u^fieE^ (484a) 



VxVxH<. = Vx (-iweEc) 
= —iLoeV X Er 



We introduce the wavenumber k as 



— ici;£:(ici;/iHc) 

wV^He (484b) 



Lo^He (485) 



Applying (560j) to (484) results in 

V X V X Ec = V(V ■ Ec) - V ■ (VEc) = u^fieE^ 
V(0) - V ■ (VEJ + w VEc = 

V^Ec + k^Ec = (486a) 
V^Hc + k^Hc = (486b) 

where the derivation of (486b) follows (486a). Thus both the electric and magnetic fields satisfy 
the vector wave equation (486a)-(486b), as well as the Maxwell equations for time-harmonic fields 
(483). In fact, any two vector fields M and N which satisfy (486) and (483) are valid electric and 
magnetic fields. In keeping with convention [Bohren and Huffman, 1983], we let M and N denote 
a candidate pair of solutions to the electro-magnetic field equations. 

We shall apply a gauge transformation to simplify the Maxwell equations for a plane wave 
scattering from a sphere. This transformation identifies a scalar generating function which satisfies 
the scalar wave equation. This generating function then determines M and N. The problem thus 
reduces to identifying the correct scalar field. Denote the scalar field by u and introduce a constant 
vector function c such that 

M = V X (nc) (487) 
= Vm X c + mV X c (488) 

where (488) follows from the vector identity for the curl of a product (560v). The divergence of 
M (487) is zero since M is the curl of a vector field (560s) 

V-M = (489) 

The curl of M (487) is 

VxM = Vx[Vx (uc)] 

= V[V ■ {uc)]-V\uc) (490) 
= V X (Vu X c + mV X c) (491) 



7 RADIATIVE EFFECTS OF AEROSOLS AND CLOUDS 



125 



where (492) comes from the vector identity for the curl of a curl (560j) and (491) arises from direct 
subsitution of (488). Neither of these expressions is particularly appealing for direct evaluation. 
We proceed from (492) armed with the vector identities for the gradient of a divergence (560n) and 
the Laplacian of a vector (560k) 

= V X (V?i X c + uV X c) 

= V[c ■ Vm + mV ■ c] - V[V ■ (uc)] - V X [V X (uc)] 
Applying (560m) to (489) 

fxm (492) 

. . . leads to the desired result 

V^M + k^M = V X [c{V\ + k^u)] (493) 

If u satisfies the scalar wave equation then the RHS of (494) vanishes and M must satisfy the 
vector wave equation. 

V^M + k^M = (494) 

As M was constructed from the curl of c (487) , so may we construct a new vector function N 
from the curl of M 

V X M 

N = — (495) 

k 

The divergence of N vanishes since N is the curl of a vector field (560s). Hence N also satisfies 
the vector wave equation 

V^N + k^N = (496) 

This implies that 

V X N = kM (497) 

We constructed M and N to behave like electric and magnetic fields, respectively. First, they 
are construced from the curls of vector fields and so their divergence vanishes (560s), as (483a)- 
(483b). Second, relations (495) and (497) show that the curls of M and N are mutually propor- 
tional, as are (483c)-(483d). Third, relations (494) and (496) show that M and N satisfy the vector 
wave equations. 

7.2 Separation of Variables 

We now seek a function u which satisfies the scalar wave equation (the Helmholtz equation) in 
spherical coordinates 

V^M + k^u = (498) 



2^ \^ V +^^7^ h^^^T^ + 2 • 2. 7^7^ + = (499) 
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The separation of variables procedure attempts to find u as the product of three independent 
functions — one for each radial and angular coordinate: 

u{r,e,(j)) = R{r)e{e)^(j)) (500) 

where R{r), 6(6'), and $(0) contain the radial, polar, and azimuthally dependent components of 
u. Substituting (500) into (499) and dividing by i?0$ we obtain [Arfken, 1985, p. 1 15] 

1 d f o(^R\ 1 d / rie\ 1 . 



flj.2 \ dr J 9r2 sin 6 d9 \ AO J ^r^ sin 9 dcj)"^ 
The differentials are full, not partial. 

7.2.1 Azimuthal Solutions 

We isolate the second order derivative of $ by multiplying (501) by sin^ 6 — this is ill-defined 
at the poles, a detail we neglect: 



- r sm 6 



2 Id/ ^dR\ 1 d f . 



"k i?r2 dr V dr J er2 sin ^ d^ \ ^ d6 ) 



(502) 



Equation (502) has the interesting property that its LHS and RHS are equal yet depend on com- 
pletely independent variables (0 and r, 9, respectively). This can be true for all values of r, 9, 
and if and only if both sides are constant — the so-called separation constant. Since the LHS is 
an azimuthal angle, we expect $(0) to be periodic, i.e., to repeat every 2n radians. Recognizing 
the LHS allows simple trigonometric solutions, we denote the separation constant by — m^. 



1 d2$ 
$d02 

d2$ 



—rn? 



d02 

= Ci,„ei"^<^ + C2,„e-"^<^ (503) 
$m(0) = Ci^rn COS m(j) + C2,m sin m0 (504) 

where Ci m and C2,m are constants of integration determined by the boundary conditions. A re- 
quirement of our physical problem, Mie scattering, is that the resulting electromagnetic fields be 
single-valued. In other words, a given azimuthal angle must be associated with only one func- 
tional value so that 

$(0 + 271) = $(0) (505) 

Requirement (505) can be met only if m is an integer. Consider, for example m = In this 
hypothetical case, $(0) is double- valued since 

cos m0 = cos I 7^ cos = cos(| +%) = — cos | (506) 



This does not satisfy requirement (505). Hence the physical requirements of the Mie problem 
require that m is integral. 
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The symmetries involved in Mie scattering make it useful to distinguish the even and odd 
components of $. We note that cos m0 and sin mcf) are symmetric and anti- symmetric, respectively, 
so that 

Ci^m COS met) = Ci^m cos[(— m)0] (507a) 
C2,m sin m(f) = -6*2,^ sin[(-m)0] (507b) 

Following Bohren and Huffman [1983], p. 85, we subscript even functions with e and odd function 
with o. 

$e,m = Ci,mCOSm0 

^o,m = C2,m sinm0 

^ra{4>) = $e,m + $o,m (508a) 

The negative sign arising from anti-symmetric functions (507b) may be subsumed into C2,m so that 
all linear all linear combinations of $e,m and $o,m are obtainable with positive integers m. Hence 
we may assume that m takes only positive (or zero) integer values. 

Since cos m(j) and sin m0 are orthonormal, any azimuthal function can be decomposed into a 
series of these basis functions 

oo 

$(0) = cosm0 + C2,m sinm0 (509) 

m=0 



7.2.2 Polar Solutions 

The polar and radial components of (502) now must satisfy 



1 d 



,dR 
dr 



1 d / . de 

^ er2sin^d^ \ d9 



m 



sin 6 



(510) 



To separate the radial and polar terms, multiply (510) by and re-group 



1 d 
^Rd^ 



- 



1 



dr J Q sin 6 dO 



sin 6 



1 d 



de 
d^ 

,d^\ 

dr / 



m 



2„2 



1 2 2 

-k r 



+ kV 



e sin 9 d9 



sin 6 



de 
d^ 



m 



sm' 



The LHS and RHS of Equation (511) must each equal a (new) separation constant for the same 
reasons given for the separation constant in Equation (502). We label this separation constant Q. 



1 



esin^d^ 



sin^ 



de\ 



+ 



m 



1 d 

sk^ede 



sin 6 



de 

'dB 



de J sin' e 
Q 



m 



sin' 9 



e 



Q 





(512) 
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The radial component is 



d 
dr 



1 d 
^Rd^ 

dr 



dr / 

+ [kV2 



Q]R 



Q 





(513) 



It is a well-known result that the separation constant Q is quantized and takes positive integral 
values defined hy Q = n(n + 1). The reasoning for this is fxm. Substituting this into (512) 
and (513) and we obtain 



1 d 



sin 6 

d_ 
dr 



de 

d^ 



+ 



n(n + 1) 



m 



2 1 



e 







dr 



r"— I + [kV^ - n{n + 1)]R 



(514) 
(515) 



Equation (514) is the associated Legendre equation and its solutions are the associated Legen- 
dre polynomials. The solutions are Legendre polynomials of degree n and order m, denoted by 
P™(cos6'). Appendix 10.2 describes many properties of Legendre polynomials including series 
expansions and recurrence formulae. 

Legendre polynomials are usually expressed in terms of /x = cos6'. Note that the rest of this 
monograph uses /i = | cos^^| to refer to the zenith angle with respect to an entire atmosphere. The 
present section is concerned with the polar angle of a single particle, and for this use only, and to 
maintain consistency with other Mie theory presentations [e.g., Bohren and Huffman, 1983; Liou, 
2002] we define /i = cos6'. 



d 

d/i 



dQ 
d/i 



+ 



n(n 



V. 



m 



2 1 



1 



e 







(516) 



The associated Legendre polynomials of degree n and order m which satisfy (516) are denoted by 
7.2.3 Radial Solutions 

We turn now to the solution of the radial equation (522), which we re-write as 



.d^i? „ dR 



— + 2r— + [W - n(n + l)]R 
dr^ dr 







(517) 



The algebraic form of (517) simplifies if we change independent variables from r to the dimen- 
sionless p 



p = kr 
r = p/k 
dr = k^^ dp 
dp = k dr 



(518a) 
(518b) 
(518c) 
(518d) 
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2 A2 



dp2 ^ Vk 



1 



+ 



2 fP\ 
kJ 



— nin + 11 



R 



dp2 



dR 



+ 2p— + [p'~n{n + l)\R = (519) 
dp 



Equation (519) is not tractable in its present form. However, a re-defined radial function will 
transform (519) into Bessel's equation. Define Z by 



R = p-^'^Z 
Z = ^pR 



(520a) 
(520b) 



and substitute it into (519) to obtain 



dp2 



+ 2p- 



dp 



dp 



-3/2 



+ 



?,p 



-5/2 



-3/2 dZ 



-3/2 dZ , ^_i/2d'^ 



P 



2 dp 2 dp 

"l/2^ + 2pV2^ + [p2 

dp 



dp2 

n{n + l)\p-^/^Z 



3p 



-1/2 



Z_pl/2^+p3/2 

4 ^ dp ^ dp2 



"^'^ p-^/'Z + 2pi/2^ + [p2 - n(n + l)]p-^'^Z 



dp 



-Z-p- + p^— -Z + 2p- + [p^-n(n + l)]Z 



p 



,d2Z 



dZ 



A^Z 



p 



+ p- 



dZ 



dp2 dp 



[p'-{n+\f]Z 









0(521) 



Equation (521) is Bessel's equation [Watson, 1958, p. 19] for half-integral values. The solutions 
are Bessel functions of half-integral order, Z = Z_^^i^(p). Appendix 10.4 describes many prop- 
erties of Bessel functions, including modified Bessel functions, asymptotic limits, and recurrence 
formulae. 

Li terms of the radial coordinate r and wavenumber k, the solutions to (517) are 



R{n] r, k) 



1 



(522) 
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The general solution to our Helmholtz equation in spherical coordinates (499), is the product (500) 
of the radial, polar, and azimuthal solutions from (522), (514), and (504), respectively 

u{n,m;r,e,(j)) = i, (kr)i'™(cos6')(Ci,mCOsm0 + Cs^^sinm^) (523) 

vkr ' 

The physical environment of Mie scattering allows us to restrict the radial solutions Z^^ i^(kr) 
to certain subsets of Bessel functions. It is advantageous to explain the rationale for doing so now 
before the nomenclature becomes even more cumbersome as it will once we expand the incident 
wave into spherical harmonics Section 10.7 and apply the boundary conditions Section 7.2.5. As 
described in Appendix 10.4, the spherical Bessel functions Zn{kr) (571b) are suitably normalized 
re-definitions of Z^^^i^(kr). The two linearly independent sets of basis function which Zn(kr) 

represents are (1) jn(kr) and t/„(kr), and (2) hn\kr) and hn\kr). 

The first physical restriction is that the scattered electromagnetic field must be finite at the 
origin. ?/„(kr) is unbounded at the origin so j„(kr) must suffice for Mie solutions interior to the 
sphere. The second physical restriction is that the scattered electromagnetic field must be finite at 
large distances from the sphere, as r — > oo. j„(kr) is unbounded as r — )• oo whereas the Hankel 
functions are well-behaved. Hence hil^ (kr) and hn^ (kr) will be used for Mie solutions exterior to 
the sphere. 

Note that (523) depends on kr rather than r. The radius-to-wavelength ratio is physically 
meaningful in the Helmholtz equation, rather than the absolute radius or wavelength. Of course 
(523) only defines the generating function that accounts for the geometry of the scattering problem. 
In practice, the dielectric constant (i.e., index of refraction) does depend on the absolute wavelenth, 
and this information propogates into the solutions for the electric and magnetic fields which are 
coupled through the Maxwell equations. Some explications of Mie theory [e.g., Liou, 2002, p. 180] 
include refractive index dependence in the solution of the Helmholtz equation rather than waiting 
until applying the boundary conditions, as we shall do. 

7.2.4 Expansion of Plane Wave into Spherical Harmonics 

Consider the scattering of a linearly polarized plane wave (482) by a sphere of radius a centered 
at the origin of the Cartesian coordinate system. In this coordinate system, an incident plane wave 
may be written as 

Ei(r,t) = Eoexp[i(ki-r-a;t)] (524) 

where the subscript i denotes "incident". We may simplify the representation of Ei(r, t) without 
loss of generality by picking our coordinate system to align with the incident plane wave. By con- 
vention we choose (1) the incident wave to propogate toward the 2;-direction, and (2) the direction 
of polarization of the incident electric vector is parallel to the x-axis. Under our assumptions, the 
incident plane waves have no components in the j or k directions because we assumed linear po- 
larization in the i direction. Hence, Ej may be represented by the simple scalar Eq. Furthermore, 
the waves are traveling in the positive ^-direction, so the wavenumber vector kj in (524) collapses 
to the scalar k, and its dot product with position simplifies ki ■ r = kz. Finally, we assume the 
incident field is constant in time so that Ei(r, t) = Ei(r). Accordingly, (524) simplifies to 

Ei(r) = Eof^^'i (525) 

fxm 
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7.2.5 Boundary Conditions 

7.2.6 Mie Theory 

The solution of Maxwell's equations (474) for the geometry of the aerosol (usually considered to 
be spherical) in the medium of interest (e.g., air or ocean water) yields the connection between 
the index of refraction n and the single scattering properties t^, zu, and g. The complete solution 
to this important problem was derived independently by physicists Ludwig Lorenz in 1890 and 
Gustav Mie in 1908. The subject is usually called Mie theory in honor of the latter, but Lorenz-Mie 
theory would be more appropriate. In-depth discussions of Mie theory are presented in van de 
Hulst [1957], Hansen and Travis [1974], and Bohren and Huffman [1983]. 

Consider the scattering of a linearly polarized plane wave (237) by a sphere of radius r centered 
at the origin of the Cartesian coordinate system. The plane wave propogates toward the z-direction. 
The direction of polarization of the incident electric vector is along the x-axis. The subscripts i, p, 
and s refer to the incident wave, particle wave (i.e., the wave interior to the particle), and scattered 
wave, respectively. 

oo 

Ei(r, t) = Y, ^n[M« - iN^ii] (526a) 

n=l 

, oo 

Hi(r, t) = -— 5^ i?n[Mili + iNili] (526b) 

n=l 

oo 

Ep(r, t) = Y, ^n[c„MllL - ic^nN« ] (526c) 

n=l 

Hp(r, t) = V Er\drMoln + '^cJ<i^oln] (526d) 

"^^P n=l 
oo 

Es(r, t) = Y^ En[iaMSn ' KM^^J (526e) 

n=l 

, oo 

H,(r, t) = ~^Y. Eniibn^'Sn + ^nM^!,] (526f) 

where 

E„ ^ '^^^o (527) 

and Eq, the amplitude of the incident electric field, is defined in (237). The superscripts (1) and 
(3) in (527) indicate the type of spherical Bessel function into which the electric and magnetic 
fields are decomposed. The complex wavenumbers within the particle, k, and the medium, km, are 
defined by, respectively 

k = 2%np/Xo (528a) 
km = 27inm/Ao (528b) 
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The vector spherical harmonic expansions of M and N are 

Moin = cos07r„(cos^^)z„(p)^ — sin 0r„(cos 6')2;„(p)0 (529a) 

Mem = (529b) 

Noi„ = (529c) 

Neln = (529d) 

where the definition of the coordinate p at which the radial Bessel functions are evaluated depends 
on the field location. Within the particle, the radial coordinate is the radial distance times the 
internal wavenumber, while the radial coordinate used to evaluate the incident and scattered fields 
(which are outside the particle) 

rk : Internal field 
rkjn : Incident and Scattered fields 

7.2.7 Resonances 

Recent studies suggest that resonant absorption of sunlight by cloud droplets may constitute a sig- 
nificant and unaccounted-for solar energy sink in the atmosphere [e.g, Nussenzveig, 2003]. Many 
studies refer to the impact of resonances on absorptance [e.g., Chylek et ah, 1978a, b; Bennett 
and Rosasco, 1978; Bohren and Huffman, 1983; Guimardes and Nussenzveig, 1994; Markel and 
Shalaev, 1999; Mitchell, 2000; Markel, 2002; Nussenzveig, 2003; Zender and Talamantes, 2006]. 
Resolving all sharp resonances requires a resolution in size parameter x = 27rr/A (r — droplet 
radius, A — incident wavelength) of about 10"'' [Chylek et ah, 1978a]. 

This section describes application of resonance enhancement by absorbing particles such as 
soot in otherwise weakly absorbing spheres such as liquid cloud droplets. Markel and Shalaev 
[1999] present an exact theory for this enhancement. Markel [2002] shows how resonances amplify 
this enhancement. 

Markel and Shalaev [1999] present a general theory of energy absorption by absorbing spheres 
in weakly-absorbing media. The theory has three main assumptions: 

1. fxm 

The absorption enhancement factor G is the ratio of absorption a by an inclusion inside a relatively 
weakly absorbing medium to the absorption 0;*^°^ of the same inclusion in a vacuum. For concrete- 
ness, we may henceforth refer to the absorbing inclusion as the soot cluster and the relatively 
weakly absorbing spherical medium as the cloud droplet. 

G[rc, A = (531) 

aW(rn, A) 

where subscripts c and n refer to the cloud droplet and the soot inclusion, respectively. Hence 
is the cloud droplet radius, is the soot inclusion radius, and the cloud droplet size parameter 
Xc = 2nrc/\. 

With the above assumptions, the enhancement factor G (531) depends only on the wavelength 
and the size of the cloud droplet, not on the size of the inclusion. 
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7.2.8 Optical Efficiencies 

Mie theory predicts the optical efficiencies of particles of a given size at a given wavelength. 
Here x stands for a, s, or e which represent the processes of absorption, scattering, and extinction, 
respectively. The optical efficiency for each of these processes is the ratio between a particle's 
effective cross-sectional area for the specified interaction (absorption, scattering, or both) with 
light and its geometric cross sectional area. 

Qa(r,A) = ^i^ (532a) 

Q,(r,A) = ^^ (532b) 

Qe(r, A) = (532c) 

Thus the optical efficiencies are dimensionless. In mie, the band-mean versions of Qa,, Qs, and Qe 
are named abs_f sh, sca_f sh, and ext_f sh, respectively. 

These optical efficiencies are not independent of one another. Two of the efficiency factors, 
usually the extinction efficiency Qc and the scattering efficiency Qs, are predicted directly by Mie 
theory. The third efficiency factor, the absorption efficiency Qa, is the residual that satisfies energy 
conservation 



Qa(r,A) = Qe{r,X)-Qs{r,X) (533) 



This relation states that extinction is the sum of absorption and scattering. 
As mentioned previously (92), Mie theory shows that 



lim Qe(r,A) = 2 (534) 

X-5-00 



7.2.9 Optical Cross Sections 

As mentioned above, the optical efficiencies Q^ are the ratios between a particle's effective cross- 
sectional area for interacting with light (i.e., absorbing or scattering photons) and its geometric 
cross-sectional area. The interaction cross-sections per particle are 

a{r, A) = nr^Qair, A) (535a) 
a{r, \) = nr^Qsir, \) (535b) 
k{r, A) = nr^Qeir, A) (535c) 

The relationship between the interaction cross-sections is exactly analogous to the relationship 
between the optical efficiencies (532), so that 

k{r,X) = a{r,X) + a{r,\) (536) 

The optical cross sections have dimensions of area per particle. 
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7.2.10 Optical Depths 

Using (85) we see that a column of depth Az [m] with a homogeneous particle concentration of 
N{r) [m^^] produces optical depths of 



7.2.11 Single Scattering Albedo 

The single scattering albedo is simply the probability that, given an interaction between the photon 
and particle, the particle will be scattered rather than absorbed. For a single particle size, this 
probability may easily be expressed in terms of the optical efficiencies (532), optical cross-sections 
(535), or optical depths (537) 

zu{r,X) = Qs(r,A)/Qe(r,A) 
= air,X)/kir,X) 
= rs(r, A)/re(r,A) 

Tc(r, A) 

In mie, the band-mean version of zu is named ss_alb_f sh. 



7.2.12 Asymmetry Parameter 



ra(r, A) = nr^Q^ir, \)N{r)Az 
Ts{r, A) = Tlr'^Q,{r,X)N{r)Az 
Te{r, A) = Kr^Qeir, X)N{r)Az 



(537a) 
(537b) 
(537c) 



In mie, the band-mean version of g is named asm_prm_f sh. 
7.2.13 Mass Absorption Coefficient 

It is often a reasonable approximation to neglect the effects of particulate scattering of radiation for 
wavelengths longer than about 5 |im, i.e., in the longwave spectral region. This approximation to 
longwave radiative transfer means that many longwave band models require only one parameter to 
account for the absorption (and emission) of radiation by particles, the mass absorption coefficient 

7.3 Effective Single Scattering Properties 

In nature, particles are not monodisperse but rather appear continuous size distributions. Therefore 
the single scattering properties of each particle size must be appropriately weighted and combined 
into the net or effective single scattering properties of the entire size distribution. These effective 
properties are the optical properties of an infinitely narrow spectral region. 

The size distribution, nn(r', z) [#m~^ rn^^], describes the rate of change of particle concentra- 
tion with particle size, and is a function of position, which we denote by z for height. 



(539) 



7 RADIATIVE EFFECTS OF AEROSOLS AND CLOUDS 



135 



The total particle concentration is the integral of the number size distribution 

POO 

No= nn{r)dr (540) 







Nq has dimensions of [#m^^], particle number per unit air volume. 

The total cross-sectional area Aq is the integral of the cross-sectional area weighted by the size 
distribution 

poo 

Ao= nr'^nn{r)dr (541) 



Jo 

Aq has dimensions of [m^ m~^], particle area per unit air volume. Since Maxwell's equations (474) 
are linear, the solutions to the equations, i.e., the optical efficiencies, are linear and additive. The 
appropriate weight for each property is the particle number distribution and a factor which depends 
on the particular property. 

For the rest of this section we assume that the single scattering properties of the aerosol are 
known or can be obtained. Section 9 describes a computer program which computes these proper- 
ties for arbitrary size distributions of spherical aerosols. 

7.3.1 Effective Efficiencies 

1 f°° 

Q^{z, A) = — / nr^Q^ir, X)n,,{r, z) dr (542a) 
^0 Jo 

1 1"^ 

Q,iz, A) = — / nr^Qsir, A)n„(r, z) dr (542b) 
^0 Jo 
1 



oo 



Q,{z,X) = — nr'Qe{r,X)n^{r,z)dr (542c) 
^0 Jo 

The effective efficiencies, like the fundamental optical efficiencies (542), are dimensionless. 
7.3.2 Effective Cross Sections 

The effective cross sections are the the fundamental optical cross sections (535) integrated over the 
size distribution. 

/•oo 

a{z, A) = AoQa = / 7tr2g^(r, A)nn(r, z) dr (543a) 
Jo 



oo 



a{z, A) = AoQs = / nr'Q,{r, X)n^{r, z) dr (543b) 



oo 



k{z,X) = AqQc= / '!lr^Qc{r\X)n^{r;z)dr (543c) 
Jo 

The dimensions of the effective cross-sections are [m^m^'^], particle area per unit air volume. 
Hence, they are also known as the volume absorption, scattering, and extinction coefficients. These 
units are usually expressed as [m^^] — we feel that [m^ m"^] is much clearer. 
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7.3.3 Effective Specific Extinction Coefficients 

As mentioned earlier, specific extinction is the extinction per unit mass of particle. Thus the 
specific extinction coefficients are the effective specific cross-sections (Section 7.3.2) divided by 
the total mass of particles Mq. For completeness, we list the explicit definition here: 

7 r \\ " ^oQa /o°°7rr2ga(r,A)nn(r,z)dr 

r. ,^ ^ AoQ, nr'^Qsir, \)n^ir, z) dr 

'^ = Wr^= J-frW{r,z)dr ^^'"'^ 

7 . ,x k AoQe J^%r^Qe{r,X)n^{r,z)dr 
Mz, ^) = JT = -vr- = roc 3 I — (544c) 

7.3.4 Effective Optical Depths 

_ /"OO 

fa(-2) = aAz = il>^w = / nr'^Qi,{r, A)nn(r, z) dr (545a) 

Jo 

/■oo 

fs{z) = aAz = ifj^w = Az I nr'^Qs{r, X)nn{r, z) dr (545b) 

Jo 

_ _ r'^ 
fc{z) = kAz = ^peW = Az / nr'^Qe{r, A)nn(r, z) dr (545c) 

Jo 

7.3.5 Effective Single Scattering Albedo 

w{z, A) = a{z, X)/k{z, X) = ( [ Tir'^Qs{r, A)nn(r, z) dr) / k{z, A) (546) 



\"'0 / / 

7.3.6 Effective Asymmetry Parameter 

The effective asymmetry parameter g 

~( \^ ^'^ '!^r'^g{r,X)Q^{r,X)n^{r,z)dr 

7.4 Mean Effective Single Scattering Properties 

Atmospheric radiative transfer models (with the possible exception of line-by-hne models) work by 
discretizing the spectral region of interest into a reasonable number of finite width spectral bands. 
These bands can be, and usually are, much wider that the spectral structure of the absorption and 
scattering features they contain. For example, many GCMs divide the solar spectrum into about 
twenty bands, and the infrared spectrum into about ten'^. The appropriate optical properties for 
such finite bands are called the mean effective single scattering properties. 



'^CAM uses nineteen solar bands and six-to-eight infrared bands. The Malkmus narrow band model SWNB2 covers 
the solar spectrum with 1690 bands. 
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Naturally, we construct the mean effective properties as a weighted mean of the effective sin- 
gle scattering properties (Section 7.3). For the single scattering properities of atmospheric size 
distributions, the appropriate spectral weight is the specific flux or radiance. For solar or infrared 
radiative transfer, weight by the fractional solar or infrared flux, respectively. These procedures 
are sometimes called Rayleigh weighting or Planck weighting, respectively If possible, use the 
ambient temperature for determining infrared weights. 

Other potential spectral weighting factors depend on the application. For example, the single 
scatter properties of the surface particle size distribution determine the surface reflectance. The 
spectral surface reflectance is a more optimal weight than the incident specific flux for determining 
the band average reflectance of snow and ice surfaces. 

We denote the spectral weighting function by W{\). Spectral flux weighting in Earth's atmos- 
phere may be done by setting 

where F^(\) is the solar spectral irradiance and Bx(\, T) is the specific blackbody radiance (41b). 
For brevity, we omit the explicit dependence of W{\) on any factor (e.g., altitude, atmospheric 
composition) except wavelength. The mean effective single scattering properties must be normal- 
ized by the integral of the spectral weighting function over the region of interest, Wq. 

Wq= j W{\)A\ (549) 

Unless the spectral weighting function is normalized, i.e., integrates to unity over the region of 
interest, the computation of (549) should be done once, outside the spectral integration loop. 

7.4.1 Mean Effective Efficiencies 

The mean effective efficiencies are the flux-weighted opacities of a size distribution of particles. 

^a(^) = -TTiT / / ^^'Qa(r, A)r2,(r, z)W{\) dr dA (550a) 



'^min ^ 



Q-^^) = -TUT- / / ^^'Qs(r, X)n^{r, z)W{X) dr dA (550b) 

AoWo JA^in Jo 

_ 1 rAmax rCO 

^e(^) = Tli7/ / T^r^Qc{r,X)n^{r,z)W{X)drdX (550c) 
^ovVo Jx- Jo 



The mean effective efficiencies, like the effective efficiencies (542) and the fundamental opti- 
cal efficiencies (532), are dimensionless. In mie, Qa., Qs, and Qc, are named abs_f sh_f f c, 
sct_f sh_f f c, and ext_f sh_f f c, respectively. 
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_ 1 /"Amax rOO 

a{z) = AoQ, = — / 7rr2g,(r, A)n„(r, z)W{X) dr dX (551a) 

yyo Jx^,„ Jo 



_ 1 /"Amax poo 

a{z) = AoQs = T7r / 7Cr'a(r, A)n,(r, z)W{\) dr d\ (551b) 

_ _ 1 /"Amax /-OO 

fc(^) = AoQe = 777/ / 7Cr2Qe(r, A)n,(r, ^)iy (A) dr dA (551c) 

^0 Jx,,,in Jo 

The mean effective cross sections, like the effective cross section (543), and the fundamental opti- 
cal cross sections (535), have dimensions of area. In mie, a, a, and k, are named abs_cf f _vlm, 
sca_cf f _vlm, and ext_cf f _vlm, respectively. 

7.4.3 Mean Effective Specific Extinction Coefficients 

. . ^ ^ ^ lt:r^r^Q^ir,>^W{r,z)WiX)drdX 

J,'-W{X)dXj-fr^n^{r,z)dr 

^ — « J,-nr^Q.{r,XW{r,z)W{X)drdX 

Wsiz) = — = = — (552b) 

j;:-W{X)dXj-fr^n^{r,z)dr 

, , , ^^.-.e :t:Io^^r^Qe{rA)n^ir,z)WiX)drdX 

ipeiz) = — = = — (552c) 

J^-W{X)dXj-fr^n^{r,z)dr 

The mean effective specific extinction coefficients, like the effective specific extinction coefficients 
(544) have dimensions of area per unit mass. In mie, and upe are named abs_cf f _mss, 

sca_cf f _mss, and ext_cf f _mss, respectively. 

7.4.4 Mean Effective Optical Depths 

/-Amax rOO 

fa = aAz = ^^w = — / nr^Q^ir, X)n^{r, z)W{X) dr dX (553a) 

VVo JA^in Jo 

f, = aAz = ijj,w= — / 7rr2Qs(r, A)nn(r, z)W{X) dr dX (553b) 



a 


AoQ^ 


Mo 


Mo 


a 


AoQs 


Mo 


Mo 


k 


AoQc 


Mo 


Mo 



Wo 



Amin J 



fe = kAz = = / ^r^Qeir, A)nn(r, z)W{X) dr dX (553c) 

Wo Jx^_^.^_^ Jo 

7.4.5 Mean Effective Single Scattering Albedo 

1 /"Amax /"CO 

w{z) = a{z)/k{z) = —J— / / nr'Q.ir, A)n,(r, z)W{X) dr dX (554) 
Wo/c(z) Jx^,^ Jo 

In mie, is named ss_alb. 
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7.4.6 Mean Effective Asymmetry Parameter 



The mean effective asymmetry parameter g 



9[z) = 




(555) 



In mie, g is named asm_prm. 

7.5 Bulk Layer Single Scattering Properties 

The real atmosphere is characterized by multiple species coexisting and interacting with the radi- 
ation field. The radiation field is determined by the combined optical properties of all the radia- 
tively active constituents. Knowing only the individual radiative properties of the constituents is 
not helpful as the radiation field is not in any sense linear, i.e., the additive result of the radiation 
fields produced by each constituent individually. Instead the individual radiative properties of all 
constituents, gases, particles, and boundary surfaces must be combined into what we shall call 
the bulk single scattering properties or layer single scattering properties. Combination of optical 
properties proceeds as in the previous section. 

7.5.1 Addition of Optical Properties 

Let denote the number of radiatively active species in a volume. If we assume that these species 
do not interact with eachother then the scattering and absorbing properties of the medium are 
additive. This may be called the independent scatterers assumption. The combination rules that 
result from this assumption were described by Cess [1985]. 

The independent scatterers assumption has many critics. Melnikova [2008] believes that enough 
interactions occur on the scale of interstitial aerosol and gases in clouds, that the independent scat- 
terers assumption significantly underestimates the effective optical depth (and thus absorption) of 
these media. Knyazikhin et al. [2002] may also question the assumption. 

7.5.2 Bulk Optical Depths 

Optical depths add linearly. With reference to (545), 



N 




(556a) 



i=l 



N 




(556b) 



i=l 



N 




(556c) 



i=l 
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7.5.3 Bulk Single Scattering Albedo 



With reference to (546), (545) and (556c), 



w = — 



1 



N 



(557) 



7.5.4 Bulk Asymmetry Parameter 



With reference to (555), (546), (545), (556c), and (557), 



9= 




(558) 



7.5.5 Diagnostics 

When the radiation field is known, it is possible to diagnose the contribution of individual elements 
to bulk properties such as heating. For instance, the contribution of given- size particles to the 
net radiative heating of a layer is the difference between the absorption and the emission of the 
particles. 



8 Global Radiative Forcing 

Figure 3 shows the balance between incident, reflected, and absorbed solar radiation observed by 
the ERBE satellite system from 1985-1989. 

Figure 4 shows the balance between terrestrial radiation emitted by the surface, trapped by the 
atmosphere, and escaping to space as observed by the ERBE satellite system from 1985-1989. 

Figure 5 shows the temperature and OLR response to ENSO in terms of HovmoUer diagrams. 

Figure 6 shows the shortwave cloud forcing (SWCF) observed by the ERBE satellite system 
from 1985-1989. 

Figure 7 shows the zonal mean shortwave cloud forcing (SWCF) observed by the ERBE satel- 
lite system from 1985-1989. 

Figure 8 shows the longwave cloud forcing (LWCF) observed by the ERBE satellite system 
from 1985-1989. 

Figure 9 shows the zonal mean longwave cloud forcing (LWCF) observed by the ERBE satellite 
system from 1985-1989. 

Figure 10 shows the balance between solar and terrestrial cloud forcing as observed by the 
ERBE satellite system from 1985-1989. 

Figure 1 1 shows the interannual variability of the cloud forcing over the Pacific in terms of 
HovmoUer diagrams of cloud forcing. 




(559) 
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Figure 3: Geographic distribution of 1985-1989 climatological mean (a) insolation cJid'^-Fsw (b) 
reflected shortwave irradiance dd'^-^sw ^^'^ absorbed shortwave radiation Ji^^F^w [Wm^^] 
from ERBE observations. 
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lRBE Surface-Emitted Longwave Flux (W m ^) 81 
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Figure 4: Geographic distribution of 1985-1989 climatological mean terrestrial radiation (a) emit- 
ted by the surface cm-^lw (b) trapped by the atmosphere Gc, and (c) escaping to space ci^^F^-ij^ in 
[Wm~^] from ERBE observations. 




135"E 180" 135'W 90'W 1 35'E 180' 1 35'W 90'W 



Figure 5: HovmoUer diagrams of (a) sea surface temperature [K] and (b) outgoing longwave 
radiation [Wm^^] over the Equatorial Pacific (averaged 10°S-10°N). Month 1 is January 1985. 
Contour intervals are 0.5 K and 10 Wm^^, respectively. 
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Figure 6: Geographic distribution of shortwave cloud forcing SWCF [Wm ^] for 1985-1989 
from ERBE observations for (a) January and (b) July. 
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Figure 7: Zonal mean shortwave cloud forcing SWCF [Wm ^] for 1985-1989 from ERBE ob- 
servations and from CCM simulations for (a) January and (b) July. 
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ERBE Longwave Cloud Forcing (W m ^) 
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Figure 8: Geographic distribution of shortwave cloud forcing LWCF [Wm ^] for 1985-1989 
from ERBE observations for (a) January and (b) July. 
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Figure 9: Zonal mean shortwave cloud forcing LWCF [Wm ^] for 1985-1989 from ERBE ob- 
servations and from CCM simulations for (a) January and (b) July. 
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Figure 10: Geographic distribution of 1985-1989 climatological mean (a) shortwave cloud forcing 
SWCF, (b) longwave cloud forcing LWCF, and (c) net cloud forcing NCF [W m^^] from ERBE. 
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Figure 11: HovmoUer diagrams of cloud forcing [Wm^^] in the Equatorial Pacific (averaged 
10°S-10°N). ERBE observations of (a) shortwave cloud forcing SWCF, and (b) longwave cloud 
forcing. Month 1 is January 1985. Contour interval is 10 Wm^^. 
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9 Implementation in NCAR models 

The discussion thus far has centered on the theoretical considerations of radiative transfer. In 
practice, these ideas must be implemented in computer codes which model, e.g., shortwave or 
longwave atmospheric fluxes and heating rates. This section describes how these ideas have been 
implemented in the NCAR narrow band models. 



10 APPENDIX 

10 Appendix 



151 



10.1 Vector Identities 

A number of vector identities are useful in the derivation of Mie theory. For arbitrary vectors a, b, 
and c; arbitrary scalars a and j3; and scalar fields u = u{x, y, z) and v = v{x, y, z), we have: 



Grad-linearity: 


V{au + (3v) 


= avu + pvv 




(560a) 


Div-hnearity: V ■ (aa + /3b) 


= aV • a + p V • b 




(SoUb) 


Curl-linearity: V x (aa + /3b) 


= aV X a + pV X b 




(5dUc) 


Dot-commutativity: a ■ b 


1 

= b ■ a 




(560d) 


Cross-(non)-commutativity: a x b 


= — b X a 




(SoUe) 


Scalar Triple Product: 


a ■ b X c 


=caxb=bcxa 




(5601) 


Cross-cross: 


a X (b X c) 


= b(a ■ c) — c(a • b) 




(560g) 




(a X b) X c 


1 / \ /I \ 

= b(a ■ c) — a(b ■ c) 




(560h) 


Curl-cross: 


V X (a X b) 


= a(V ■ b) - b(V ■ a) + (b • 


■ V)a-( 


> ■ V)b (560i) 


Curl-curl: 


V X (V X a) 


= V(V ■ a) - V ■ (Va) 










= V(V ■ a) - V^a 




(560j) 


Vector Laplacian: V^a 


= (V ■ V)a = V ■ (Va) 




(560k) 






= V(V ■ a) - V X (V X a) 




(5601) 


Grad-div: 


V(V ■ a) 


= V^a + V X (V X a) 




(560m) 






7^ (not equal) V a 






Grad-dot: 


V(a-b) 


= a X (V X b) + b X (V X , 


a) + (b ■ 


V)a+ (a- V)b 






— [SL ■ v )D -\- [D ■ vja + ax 


(V X b 


j + D X \y X aj 










(560n) 


Del-vector: Va 


= VifliXj = Dyadic 




(560o) 


Div-grad: 


V ■ (Vn) 


= V\ 




(560p) 


Div-del: 


(V ■ V)a 


= V^a 




(560q) 


Dot-cross: 


a ■ (b X c) 


= b • (c X a) = c ■ (a X b) 




(560r) 


Div-curl: 


V ■ V X a 


= 




(560s) 


Curl-grad: 


V X Vm 


= 




(560t) 


Div-product: 


V ■ (wa) 


= mV ■ a + a ■ Vm 




(560u) 


Curl-product: 


V X (wa) 


= Vm X a + uV X a 




(560v) 


Dot-del: 


(a- V)a 


= ^V|ap + (V X a) X a 




(560w) 



The divergence of the curl vanishes (560s). The curl of the gradient vanishes (560t). 

The Laplacian V^ is shorthand for (V ■ V), and is usually read as "del-squared", though it is 
not the "square" of the del operator. The scalar Laplacian is the divergence of the gradient (560p). 
However, the vector Laplacian is not the gradient of the divergence (560m). The vector Laplacian 
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is defined by (??) which is obtained by simply re-arranging (560j). The vector Laplacian (??) 
is the gradient of the divergence (560m) minus the curl of the curl (560j). 
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10.2 Legendre Polynomials 

The associated Legendre equation of degree n and order m is 



153 



d 

dyU 









+ 



n(n + 1) 



m 



e 







(561) 



The solutions to (561) are associated Legendre polynomials, also called spherical harmonics. The 
associated Legendre polynomial of degree n, order m, and argument /i is denoted by P^{jj,). The 
P™(/i) are orthogonal functions 



[n + m)! 



{n — m)! 



(562) 



When m = 0, P^{jJ,) reduces to the Legendre polynomial of order n, Pn{l^) whose orthogonality 
properties follow from (562) and are shown in (141). Applying suitable normalization to P^{^) 
yields the normalized associated Legendre polynomials A™(/i) 



An(/^) 



2n + 1 (n — m)! 



[n + 171)1 



(563) 



The A™(/i) are orthonormal and are the polar components of the spherical harmonics Y^{9, cj)) 
(566). 

10.3 Spherical Harmonics 

Spherical harmonics are the product of the normalized associated Legendre polynomials (563) and 
the orthonormal exponentials (564) 



\m<l> 



2n 



(564) 



The spherical harmonics Y"' 



of degree n and order m are 



n 



= (-V. 



) ^P^(cos6')e' 



(565) 



The Condon-Shortley phase factor (— 1)™ is a convention which causes the sign alternation for the 
positive m spherical harmonics [Arflcen, 1985, p. 682]. 
The spherical harmonics (566) are normalized such that 



271 rd=-K 



(566) 



=0 Je=o 



where * denotes the complex conjugate. 
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10.4 Bessel Functions 

The definitive reference for properties of Bessel functions is Watson [1958]. Other useful refer- 
ences are Abramowitz and Stegun [1964], p. 374 and Arjken [1985], p. 573. When the Helmholtz 
equation is solved by separation of variables in cyhndrical or spherical coordinates (521), the radial 
component satisfies Bessel's equation [Watson, 1958, p. 19], 

z^^ + z— + [z^ - u^]w = (567) 
az^ az 

Bessel functions are solutions to (567). The order of the functions is indicated by the parameter 
u (not to be confused with frequency). The functions take a single argument z. Both u and z may 
be complex, e.g., z = x + iy, although our particular applications depend only on real- valued u. 
Indeed, we require only integer and half-integer v. 

Due to their power and range of applicability, a menagerie of inter-related solutions to (567) 
exist. Fortunately, the terminology to identify and describe them is standard. The standard Bessel 
functions of the first kind, second kind, and third kind, are denoted by Ju{z), Yi,{z), and hI^\z) 
and Hy {z), respectively. Bessel functions of the second kind are also called Neumann functions 
and notated with is N^{z) instead of Y,,{z). Mathematical physics texts often use the Neumann 
function notation N^{z) to avoid confusion with spherical harmonics; whereas pure mathematical 
texts normally use the Y^(z) notation. Bessel functions of the third kind are also called Hankel 
functions. All of these satisfy (567) and so are properly called Bessel functions. Specialized ver- 
sions of Bessel functions derive from the three basic kinds. These specializations include modified 
Bessel functions, and spherical Bessel functions. 

General solutions to the Bessel equation (567) are formed from any two linearly independent 
Bessel functions. For non-integer u, Ju{z) and J-^{z) are linearly independent solutions and so 
form a complete basis. However, for integral u 

J_n{z) = i-lTJniz) Integern (568) 

and these solutions are not linearly independent. 

Bessel functions of the second kind Y^^z) are always linearly independent of Jy{z). They may 
be defined as 

Y^{z) = N^^,^ = Mz)cos{un)-J..{z) ^^^^^ 

sin un 

where the RHS limiting value used for integral u. Bessel functions of the first and second kind, 
Juiz) and Y^{z) together form a complete basis. Series expansion (not shown) reveals that Yn{z) 
diverges logarithmically as z ^ 0. Boundary conditions which require finite solutions at the origin 
(z = 0) therefore automatically exclude Yn{z). Conversely, Yn{z) may appear in any solution 
which does not require finite values at the origin. 

Bessel functions of the third kind (Hankel functions), H^y\z) and hI^\z), are another set of 
linearly independent solutions to Bessel's equation. The relation between these complete bases is 

Hl'\z) = J,{z) + iY,{z) (570a) 
Hi^\z) = Mz)-iY,{z) (570b) 



10 APPENDIX 155 

10.4.1 Spherical Bessel Functions 

The spherical Bessel functions of the first and second kind are 

U^) = \lf/.^^k^^) (571a) 

l/^W = /^n+i;^W (571b) 

In analogy to (570b) the spherical Bessel functions of the third kind, or spherical Hankel func- 
tions are 

h^i^\z)=Uz) + iy,{z) (572a) 
h^;^\z)=Uz)-iy,{z) (572b) 

10.4.2 Recurrence Relations 

Bessel functions (and thus spherical Bessel functions) satisfy the recurrence relations 

Z,^i{z) + Z,+i{z) = '^Z^iz) (573a) 
Z,_i(^) - Z,+,iz) = 2Zliz) (573b) 
Zliz) = Z,^i{z) - -^Z,{z) (573c) 

Zl{z) = -Z,+i{z) + -^Z,{z) (573d) 
Seed values for recurrence relations (573d) may be obtained from 



Jo{z) = -Ji{z) (574a) 
Fo'(^) = -Yi{z) (574b) 

10.4.3 Power Series Representation 

The power series for regular Bessel functions of the first kind is 

fc=0 ^ ^ 

Modified Bessel functions of the first kind are defined in terms ordinary Bessel functions as 

Uz) = [-".Miz) (576) 

The power series for the modified Bessel function of the first kind, luiz) differs from (575) due to 
the absence of the negative sign 

fc=0 ^ ^ 
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Of particular interest to radiative transfer are Io{z) and Ii{z) which appear in the solution to the 
mean beam absorptance (397) and equivalent width of an isolated Lorentz line. 

oo /I \2k 



..0 *! (* + 1)' 



For z <^ 1, (578) becomes 



/o(2 < 1) ^ 1 (579a) 
h{z < 1) ^ z/2 (579b) 

10.4.4 Asymptotic Values 

The asymptotic behavior of functions can often be ascertained from their power series representa- 
tion. Bessel functions obey the following limits as 2; — > 0: 

Z,{z) ^{\zY/T{v+l) (z/^ -1,-2,-3...) (580a) 
Yq{z) {2/n)hiz (580b) 

^ (580C) 
10.5 Gaussian Quadrature 

The accuracy of automatic integration of analytic formulae by computers rests on three criteria: 

1 . Smoothness of function 

2. Number of quadrature points 

3. Placement of quadrature points 

For many problems the functional form of the integrand is given. We are usually free to change, 
however, the number and placement of quadrature points. For many applications, the optimal 
number and location of quadrature points are related through Gaussian quadrature. 
Gaussian quadrature points Xk have the property that 

/b " 
/(x)^i;(x)dx^5^Afe/(xfc) (581) 
k=i 

is exact if /(x) is any polynomial of degree < 2ri — 1. The values of the Xk are determined by 
the form of the weighting function w{x) and the interval boundaries [a, h]. In fact the Xk are the 
roots of the orthogonal polynomial with weighting function w{x) over the interval [a, h] [see, e.g., 
Arfken, 1985, p. 969]. A common case in radiative transfer has w{x) = 1, [a, b] = [—1, 1] 

/I 
/(x)dx^5^A,/(xfc) (582) 
1 k=i 
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For this case, the Xk are the n roots of the Legendre polynomial of degree n, Pn{x), and the 
Gaussian weights Ak are defined by the derivative of Pn{x) 

Xk = xe [0, 1] I Pn{xk) = (583) 

Radiative transfer problems often use Gaussian weights in the angular integration. Depending 
on whether we are concerned with net or hemispherical quantities the angular integration will have 
different intervals 

/I " 
/(xV(x)dx^^Afc/(x,) (585) 
"1 k=i 

„1 n 

/ /(i/)«;(i/)di/^ V5fc/(y,) (586) 
•^0 k=l 

Equation (585) is known as full-range Gaussian quadrature and (586) as half-range Gaussian 
quadrature. 

The half-range points and weights are intimately related to the full-range points and weights. 
The interval y E [0, 1] is a simple linear transformation of x G [—1,1] 

y = ^ (587) 
x = 2y-l (588) 

The half -range points and weights may be obtained by applying the forward transformation (587) 
to (585). Thus 

Xi. + 1 

Vk = (589) 

The weights A^ depend only on n (fxm: I think) and are independent of f(x) and w{x). For 
Legendre polynomials, the A^ and sum to the size of the interval 



n 



Y,Ak = 2 (590) 

k=l 
n 

J2Bk = l (591) 



k=l 

A common technique is to apply half of the full range quadrature points and weights to half of 
the full range interval, [0, 1]. In this case the quadrature is no longer Gaussian (i.e., not exact for 
polynomials of degree < 2n), but the quadrature points are often quite optimal. 

Table 9 lists the first five Legendre polynomials, their roots, the associated full Gaussian angles, 
and their Gaussian weights. Abramowitz andStegun [1964], p. 917, list Gaussian quadrature points 
and weights for n <32. 

Table 10 summarizes the same information needed for half -range Gaussian quadrature. 



'^When w{x) ^ 1, the Gaussian quadrature points are the roots of other orthogonal polynomials, e.g., Chebyshev, 
Laguerre, Hermite. 
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Table 9: Full-range Gaussian quadrature' 



n, h 


Polynomial 


Root 




Decimal 


Angle 


Weight 

Ak 


1,1 









0.0 


90.0° 


2.0 


2,1 


i(3/.2 - 1) 


1/v^ 




0.577350 


54.7356° 


1.0 


3,1 


i(5/i3 -3/i) 







0.0 


90.0° 


0.888889 


3,2 




Ti 




0.774597 


39.2321° 


0.555556 


4,1 


|(35;U^ - 30/^2 + 3) 


( 15-2^30^ 
V 35 , 


f 


0.339981 


70.1243° 


0.652145 


4,2 




( 15+2v^^ 
I 35 , 


r 


0.861136 


30.5556° 


0.347855 



"The /ife and Ak are only supplied for k < The /ifc are antisymmetric, and the Ak are symmetric, about 

k = n/2. Thus = -fik and A^-k = Ak- 



Table 10: Half- range Gaussian quadrature' 



n, k 


Polynomial 


Root 


Decimal 


Angle 


Weight 




^n(/i) 




l^k 


dk 


Ak 


1,1 


/i 





0.0 


90.0° 


2.0 


2,1 


l(3/i2-l) 


+ 1/^3 


+0.57735 


+54.7356° 


1.0 


3,1 


1 (5/^3 -3/.) 










4,1 


l(35/i^-30/i2 + 3) 











"The symmetries differ from Table 9. The polynomials, roots and weights need to be transformed from their 
full-range counterparts. 



10 APPENDIX 



159 



10.6 Gauss-Lobatto Quadrature 

Lobatto quadrature, is a type of Gaussian quadrature that prescribes abscissae at the integration 
endpoints. The Lobatto quadrature rule of order L is 

r f{x)dx = (592) 
"^-^ 1=1 

l=L 

= HJi-l) + Y,Hif{xi) + HLfil) (593) 
1=1 

Whereas Gaussian quadrature is exact for polynomials of degrees < 2L — 1, Lobatto quadrature 
fixes two quadrature points and so is exact only for polynomials of degrees < 2L — 3. Nevertheless, 
the placement of abscissae can give Lobatto quadrature a higher effective accuracy for certain 
classes of functions such as asymmetric phase functions [Wiscombe, 1977]. 

The free (non-prescribed) Lobatto abscissae for / = 2, . . . , L — 1 are the L — 2 zeros of the 
derivative of the Legendre polynomial of degree L — 1 

P'l-M) = (594) 

(By comparison, the abscissae in Gauss-Legendre quadrature are Legendre polynomial roots). 
The interior Lobatto weights Hi satisfy 

= L(L-l)[Vi(x.)]^ ^'''^ 



while the endpoint weights are 



Michels [1963] describes numerical approaches to determining Lobatto quadrature abscissae 
endpoints and weights. This procedure converges remarkably quickly and smoothly. However, the 
traditional quadrature variable u = cos{9) and range [—1, 1] produce abscissae that a forward- 
clustered in cos(^^) and thus angles that are relatively evenly distributed (rather than forward- 
clustered) in 9. 

Wiscombe [1977] puts forward a solution to this problem. He changes the quadrature variable 
and domain to 9 and [0, n], respectively. Implementing this procedure is relatively straightforward. 
The Lobatto weights Hi must be normalized on the new integration domain 



while the endpoint weights are 
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10.7 Exponential Integrals 

The exponential integrals En{x) are defined by 



oo 

—n„—tx 



En{x) = J t""e-*^dt (599) 

and by the behavior of at x = 

Ei{0) = oo 
EniO) = l/(n-l) 

A definition equivalent to (599) is obtained by making the change of variable w = t^^. This maps 
t G [l,oo) into w G [1,0] 

w = 
t = w'^ 
dt = —w~'^ dw 

dw = -r^dt (600) 



Substituting (600) into (599) we obtain 



E^{x) = I (w;-^)-"e-^/'"(-u;-^)dM; 



1 

w^w-'e""/" dw 

1 

^n-2^-x/w^^ (601) 



Exponential integrals satisfy two important recurrence relations 

-E'^{x) = -En-i{x) (602a) 
nE(„+i)(x) = e"" - xEn{x) (602b) 

The angular integral of many radiometric quantities may often be expressed in terms of the first 
few exponential integrals. 

E„(r) = I /x"-2e-"/^ d/i (603a) 







Ei(r) = / /i-^e-"/'^ d/i (603b) 







E2{t) = / e-"/^d/i (603c) 







E3(r) = / /ie-"/^d;U (603d) 
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plane albedo, 51 

plane-parallel, 5 

planetary albedo, 2, 52 

planetary radiative equilibrium, 1, 16 

Poisson distribution, 74, 75 

power series, 21, 21, 155 

pressure broadening, 70, 74 

pressure-broadened line width, 74 

pressure-shift, 95 

quantum efficiency, 9 
quantum yield, 9 

R-branch, 87 
radiance, 4 

radiative equilibrium, 35 

radiative equilibrium temperature profile, 35 

radiative flux divergence, 36 

radius of curvature, 120 

random band, 116 

Rayleigh limit, 62 

Rayleigh weighting, 137 

Rayleigh- Jeans limit, 17 

reciprocity relations, 46 

reduced mass, 82 

reflectance, 13, 48 
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reflection, 48 
residue theorem, 23 
resonances, 132 
resonant absorption, 132 
Riemann zeta function, 21, 21 
rigid rotator, 82 
rotational constant, 83 
rotational energy, 82 
rotational partition function, 88 
rotational quantum number, 83 

scaled absorber path, 117 

scaled pressure, 117 

scaled temperature, 117 

scaling approximation, 117 

scattered wave, 131 

scattering, 37 

scattering angle, 37, 64 

scattering efficiency, 133 

selection rule, 84 

selection rules, 86 

self-broadening, 94 

separation constant, 126 

separation of variables, 126 

shaded pyranometers, 69 

shadowband radiometers, 69 

single scattering albedo, 121 

single scattering properties, 121 

single- scattering approximation, 67 

single-scattering source function, 43, 44, 67 

size parameter, 62 

size parameters, 63 

solar constant, 2, 24 

sohd angle, 3 

source function, 25 

specific extinction, 136 

specific heat capacity at constant pressure, 36 

specific intensity, 4 

spectral, 15 

spectral absorption cross section, 70 

spectroscopy, 81 

specular, 41 

speed of light, 3 

spherical albedo, 2, 52 

spherical Bessel function, 131 
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spherical Hankel functions, 155 
spherical hannonics, 153 
spherical polar coordinates, 6, 38 
spontaneous decay, 72 
spontaneous emission, 72, 91 
spring constant, 85 
statistical weight, 88 
Stefan-Boltzmann constant, 23 
Stefan-Boltzmann Law, 2, 23, 34 
stimulated absorption, 91 
stimulated emission, 92 
stratified atmosphere, 5, 28 
strong line limit, 100, 101 
strong-line limit, 101 
surface reflectance, 29, 42, 49, 137 
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Tesla, 122 

the mass extinction coefficient, 25 
the number extinction coefficient, 25 
the volume extinction coefficient, 25 
thermal emission, 15 
thermal speed, 75 
thermodynamic equilibrium, 16 
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tranlational kinetic energy, 82 

transmission, 48 

transmittance, 48, 48 

two level atom, 91 

two stream equations, 44 

two-stream approximation, 34, 54 

upwelling, 19 

vector identities, 151 
vector spherical harmonic, 132 
vector wave equation, 124, 125 
vibration-rotation bands, 87, 88 
vibrational constant, 86 
vibrational energy, 82 
vibrational ground state, 86 
vibrational partition function, 89 
vibrational quantum number, 86 
Voigt line shape, 70 
Voigt profile, 81 
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volume extinction coefficient, 27 

warming rate, 35 
wavelength, 3 
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weak-line limit, 99 
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